!/*****************************************************************************/
! *
! *  Elmer, A Finite Element Software for Multiphysical Problems
! *
! *  Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
! *
! *  This library is free software; you can redistribute it and/or
! *  modify it under the terms of the GNU Lesser General Public
! *  License as published by the Free Software Foundation; either
! *  version 2.1 of the License, or (at your option) any later version.
! *
! *  This library is distributed in the hope that it will be useful,
! *  but WITHOUT ANY WARRANTY; without even the implied warranty of
! *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
! *  Lesser General Public License for more details.
! *
! *  You should have received a copy of the GNU Lesser General Public
! *  License along with this library (in file ../LGPL-2.1); if not, write
! *  to the Free Software Foundation, Inc., 51 Franklin Street,
! *  Fifth Floor, Boston, MA  02110-1301  USA
! *
! *****************************************************************************/
!
!/******************************************************************************
! *
! *  Authors: Juha Ruokolainen
! *  Email:   Juha.Ruokolainen@csc.fi
! *  Web:     http://www.csc.fi/elmer
! *  Address: CSC - IT Center for Science Ltd.
! *           Keilaranta 14
! *           02101 Espoo, Finland
! *
! *  Original Date: 01 Oct 1996
! *
! ****************************************************************************/

!> \ingroup ElmerLib
!> \{

!-----------------------------------------------------------------------------
!>  Module defining Gauss integration points and
!>  containing various integration routines.
!-----------------------------------------------------------------------------
MODULE Integration
   USE Types
   USE LoadMod
   USE LinearAlgebra

   IMPLICIT NONE

   INTEGER, PARAMETER, PRIVATE :: MAXN = 13, MAXNPAD = 16 ! Padded to 64-byte alignment
   INTEGER, PARAMETER, PRIVATE :: MAX_INTEGRATION_POINTS = MAXN**3

   LOGICAL, PRIVATE :: GInit = .FALSE.
   !$OMP THREADPRIVATE(GInit)

!------------------------------------------------------------------------------
   TYPE GaussIntegrationPoints_t
      INTEGER :: N
      REAL(KIND=dp), POINTER CONTIG :: u(:),v(:),w(:),s(:)
!DIR$ ATTRIBUTES ALIGN:64 :: u, v, w, s
   END TYPE GaussIntegrationPoints_t

   TYPE(GaussIntegrationPoints_t), TARGET, PRIVATE, SAVE :: IntegStuff
   !$OMP THREADPRIVATE(IntegStuff)
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Storage for 1d Gauss points, and weights. The values are computed on the
! fly (see ComputeGaussPoints1D below). These values are used for quads and
! bricks as well. To avoid NUMA issues, Points and Weights are private for each
! thread.
!------------------------------------------------------------------------------
   REAL(KIND=dp), PRIVATE, SAVE :: Points(MAXNPAD,MAXN),Weights(MAXNPAD,MAXN)
   !$OMP THREADPRIVATE(Points, Weights)
!DIR$ ATTRIBUTES ALIGN:64::Points, Weights
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Triangle - 1 point rule; exact integration of x^py^q, p+q<=1
!------------------------------------------------------------------------------
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: UTriangle1P=(/ 0.3333333333333333D0 /)
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: VTriangle1P=(/ 0.3333333333333333D0 /)
   REAL(KIND=dp),DIMENSION(1),PRIVATE :: STriangle1P=(/ 1.0000000000000000D0 /)

!------------------------------------------------------------------------------
! Triangle - 3 point rule; exact integration of x^py^q, p+q<=2
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(3), PRIVATE :: UTriangle3P = &
    (/ 0.16666666666666667D0, 0.66666666666666667D0, 0.16666666666666667D0 /)

   REAL(KIND=dp), DIMENSION(3), PRIVATE :: VTriangle3P = &
    (/ 0.16666666666666667D0, 0.16666666666666667D0, 0.66666666666666667D0 /)

   REAL(KIND=dp), DIMENSION(3), PRIVATE :: STriangle3P = &
    (/ 0.33333333333333333D0, 0.33333333333333333D0, 0.33333333333333333D0 /)

!------------------------------------------------------------------------------
! Triangle - 4 point rule; exact integration of x^py^q, p+q<=3
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTriangle4P = &
         (/  0.33333333333333333D0, 0.2000000000000000D0, &
             0.60000000000000000D0, 0.20000000000000000D0 /)

   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTriangle4P = &
         (/  0.33333333333333333D0, 0.2000000000000000D0, &
             0.20000000000000000D0, 0.60000000000000000D0 /)

   REAL(KIND=dp), DIMENSION(4), PRIVATE :: STriangle4P = &
         (/ -0.56250000000000000D0, 0.52083333333333333D0, &
             0.52083333333333333D0, 0.52083333333333333D0 /)

!------------------------------------------------------------------------------
! Triangle - 6 point rule; exact integration of x^py^q, p+q<=4
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: UTriangle6P = &
       (/ 0.091576213509771D0, 0.816847572980459D0, 0.091576213509771D0, &
          0.445948490915965D0, 0.108103018168070D0, 0.445948490915965D0 /)

   REAL(KIND=dp), DIMENSION(6), PRIVATE :: VTriangle6P = &
        (/ 0.091576213509771D0, 0.091576213509771D0, 0.816847572980459D0, &
           0.445948490915965D0, 0.445948490915965D0, 0.108103018168070D0 /)

   REAL(KIND=dp), DIMENSION(6), PRIVATE :: STriangle6P = &
        (/ 0.109951743655322D0, 0.109951743655322D0, 0.109951743655322D0, &
           0.223381589678011D0, 0.223381589678011D0, 0.223381589678011D0 /)

!------------------------------------------------------------------------------
! Triangle - 7 point rule; exact integration of x^py^q, p+q<=5
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: UTriangle7P = &
        (/ 0.333333333333333D0, 0.101286507323456D0, 0.797426985353087D0, &
           0.101286507323456D0, 0.470142064105115D0, 0.059715871789770D0, &
           0.470142064105115D0 /)

   REAL(KIND=dp), DIMENSION(7), PRIVATE :: VTriangle7P = &
        (/ 0.333333333333333D0, 0.101286507323456D0, 0.101286507323456D0, &
           0.797426985353087D0, 0.470142064105115D0, 0.470142064105115D0, &
           0.059715871789770D0 /)

   REAL(KIND=dp), DIMENSION(7), PRIVATE :: STriangle7P = &
        (/ 0.225000000000000D0, 0.125939180544827D0, 0.125939180544827D0, &
           0.125939180544827D0, 0.132394152788506D0, 0.132394152788506D0, &
           0.132394152788506D0 /)

!------------------------------------------------------------------------------
! Triangle - 11 point rule; exact integration of x^py^q, p+q<=6
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTriangle11P = &
    (/ 0.3019427231413448D-01, 0.5298143569082113D-01, &
       0.4972454892773975D-01, 0.7697772693248785D-01, &
       0.7008117469890058D+00, 0.5597774797709894D+00, &
       0.5428972301980696D+00, 0.3437947421925572D+00, &
       0.2356669356664465D+00, 0.8672623210691472D+00, &
       0.2151020995173866D+00 /)

   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTriangle11P = &
    (/ 0.2559891985673773D+00, 0.1748087863744473D-01, &
       0.6330812033358987D+00, 0.8588528075577063D+00, &
       0.2708520519075563D+00, 0.1870602768957014D-01, &
       0.2027008533579804D+00, 0.5718576583152437D+00, &
       0.1000777578531811D+00, 0.4654861310422605D-01, &
       0.3929681357810497D+00 /)

   REAL(KIND=dp), DIMENSION(11), PRIVATE :: STriangle11P = &
    (/ 0.3375321205342688D-01, 0.1148426034648707D-01, &
       0.4197958777582435D-01, 0.3098130358202468D-01, &
       0.2925899761167147D-01, 0.2778515729102349D-01, &
       0.8323049608963519D-01, 0.6825761580824108D-01, &
       0.6357334991651026D-01, 0.2649352562792455D-01, &
       0.8320249389723097D-01 /)

!------------------------------------------------------------------------------
! Triangle - 12 point rule; exact integration of x^py^q, p+q<=7
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: UTriangle12P = &
    (/ 0.6232720494911090D+00, 0.3215024938520235D+00, &
       0.5522545665686063D-01, 0.2777161669760336D+00, &
       0.5158423343535236D+00, 0.2064414986704435D+00, &
       0.3432430294535058D-01, 0.3047265008682535D+00, &
       0.6609491961864082D+00, 0.6238226509441210D-01, &
       0.8700998678316921D+00, 0.6751786707389329D-01 /)

   REAL(KIND=dp), DIMENSION(12), PRIVATE :: VTriangle12P = &
    (/ 0.3215024938520235D+00, 0.5522545665686063D-01, &
       0.6232720494911090D+00, 0.5158423343535236D+00, &
       0.2064414986704435D+00, 0.2777161669760336D+00, &
       0.3047265008682535D+00, 0.6609491961864082D+00, &
       0.3432430294535058D-01, 0.8700998678316921D+00, &
       0.6751786707389329D-01, 0.6238226509441210D-01 /)

   REAL(KIND=dp), DIMENSION(12), PRIVATE :: STriangle12P = &
    (/ 0.4388140871440586D-01, 0.4388140871440586D-01, &
       0.4388140871440587D-01, 0.6749318700971417D-01, &
       0.6749318700971417D-01, 0.6749318700971417D-01, &
       0.2877504278510970D-01, 0.2877504278510970D-01, &
       0.2877504278510969D-01, 0.2651702815743698D-01, &
       0.2651702815743698D-01, 0.2651702815743698D-01 /)

!------------------------------------------------------------------------------
! Triangle - 17 point rule; exact integration of x^py^q, p+q<=8
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(17), PRIVATE :: UTriangle17P = &
    (/ 0.2292423642627924D+00, 0.4951220175479885D-01, &
       0.3655948407066446D+00, 0.4364350639589269D+00, &
       0.1596405673569602D+00, 0.9336507149305228D+00, &
       0.5219569066777245D+00, 0.7110782758797098D+00, &
       0.5288509041694864D+00, 0.1396967677642513D-01, &
       0.4205421906708996D-01, 0.4651359156686354D-01, &
       0.1975981349257204D+00, 0.7836841874017514D+00, &
       0.4232808751402256D-01, 0.4557097415216423D+00, &
       0.2358934246935281D+00 /)

   REAL(KIND=dp), DIMENSION(17), PRIVATE :: VTriangle17P = &
    (/ 0.5117407211006358D+00, 0.7589103637479163D+00, &
       0.1529647481767193D+00, 0.3151398735074337D-01, &
       0.5117868393288316D-01, 0.1964516824106966D-01, &
       0.2347490459725670D+00, 0.4908682577187765D-01, &
       0.4382237537321878D+00, 0.3300210677033395D-01, &
       0.2088758614636060D+00, 0.9208929246654702D+00, &
       0.2742740954674795D+00, 0.1654585179097472D+00, &
       0.4930011699833554D+00, 0.4080804967846944D+00, &
       0.7127872162741824D+00 /)

   REAL(KIND=dp), DIMENSION(17), PRIVATE :: STriangle17P = &
    (/ 0.5956595662857148D-01, 0.2813390230006461D-01, &
       0.3500735477096827D-01, 0.2438077450393263D-01, &
       0.2843374448051010D-01, 0.7822856634218779D-02, &
       0.5179111341004783D-01, 0.3134229539096806D-01, &
       0.2454951584925144D-01, 0.5371382557647114D-02, &
       0.2571565514768072D-01, 0.1045933340802507D-01, &
       0.4937780841212319D-01, 0.2824772362317584D-01, &
       0.3218881684015661D-01, 0.2522089247693226D-01, &
       0.3239087356572598D-01 /)

!------------------------------------------------------------------------------
! Triangle - 20 point rule; exact integration of x^py^q, p+q<=9
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: UTriangle20P = &
   (/ 0.2469118866487856D-01, 0.3348782965514246D-00, &
      0.4162560937597861D-00, 0.1832492889417431D-00, &
      0.2183952668281443D-00, 0.4523362527443628D-01, &
      0.4872975112073226D-00, 0.7470127381316580D-00, &
      0.7390287107520658D-00, 0.3452260444515281D-01, &
      0.4946745467572288D-00, 0.3747439678780460D-01, &
      0.2257524791528391D-00, 0.9107964437563798D-00, &
      0.4254445629399445D-00, 0.1332215072275240D-00, &
      0.5002480151788234D-00, 0.4411517722238269D-01, &
      0.1858526744057914D-00, 0.6300024376672695D-00 /)

   REAL(KIND=dp), DIMENSION(20), PRIVATE :: VTriangle20P = &
   (/ 0.4783451248176442D-00, 0.3373844236168988D-00, &
      0.1244378463254732D-00, 0.6365569723648120D-00, &
      0.3899759363237886D-01, 0.9093437140096456D-00, &
      0.1968266037596590D-01, 0.2191311129347709D-00, &
      0.3833588560240875D-01, 0.7389795063475102D-00, &
      0.4800989285800525D-00, 0.2175137165318176D-00, &
      0.7404716820879975D-00, 0.4413531509926682D-01, &
      0.4431292142978816D-00, 0.4440953593652837D-00, &
      0.1430831401051367D-00, 0.4392970158411878D-01, &
      0.1973209364545017D-00, 0.1979381059170009D-00 /)

   REAL(KIND=dp), DIMENSION(20), PRIVATE :: STriangle20P = &
   (/ 0.1776913091122958D-01, 0.4667544936904065D-01, &
      0.2965283331432967D-01, 0.3880447634997608D-01, &
      0.2251511457011248D-01, 0.1314162394636178D-01, &
      0.1560341736610505D-01, 0.1967065434689744D-01, &
      0.2247962849501080D-01, 0.2087108394969067D-01, &
      0.1787661200700672D-01, 0.2147695865607915D-01, &
      0.2040998247303970D-01, 0.1270342300533680D-01, &
      0.3688713099356314D-01, 0.3813199811535777D-01, &
      0.1508642325812160D-01, 0.1238422287692121D-01, &
      0.3995072336992735D-01, 0.3790911262589247D-01 /)

!------------------------------------------------------------------------------
! Tetrahedron - 1 point rule; exact integration of x^py^qz^r, p+q+r<=1
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: UTetra1P = (/ 0.25D0 /)
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: VTetra1P = (/ 0.25D0 /)
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: WTetra1P = (/ 0.25D0 /)
   REAL(KIND=dp), DIMENSION(1), PRIVATE :: STetra1P = (/ 1.00D0 /)

!------------------------------------------------------------------------------
! Tetrahedron - 4 point rule; exact integration of x^py^qz^r, p+q+r<=2
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UTetra4P = &
    (/ 0.1757281246520584D0, 0.2445310270213291D0, &
       0.5556470949048655D0, 0.0240937534217468D0 /)

   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VTetra4P = &
    (/ 0.5656137776620919D0, 0.0501800797762026D0, &
       0.1487681308666864D0, 0.2354380116950194D0 /)

   REAL(KIND=dp), DIMENSION(4), PRIVATE :: WTetra4P = &
    (/ 0.2180665126782654D0, 0.5635595064952189D0, &
       0.0350112499848832D0, 0.1833627308416330D0 /)

   REAL(KIND=dp), DIMENSION(4), PRIVATE :: STetra4P = &
    (/ 0.2500000000000000D0, 0.2500000000000000D0, &
       0.2500000000000000D0, 0.2500000000000000D0 /)
!------------------------------------------------------------------------------
! Tetrahedron - 5 point rule; exact integration of x^py^qz^r, p+q+r<=3
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: UTetra5P =  &
    (/ 0.25000000000000000D0, 0.50000000000000000D0, &
       0.16666666666666667D0, 0.16666666666666667D0, &
       0.16666666666666667D0 /)

   REAL(KIND=dp), DIMENSION(5), PRIVATE :: VTetra5P =  &
    (/ 0.25000000000000000D0, 0.16666666666666667D0, &
       0.50000000000000000D0, 0.16666666666666667D0, &
       0.16666666666666667D0 /)

   REAL(KIND=dp), DIMENSION(5), PRIVATE :: WTetra5P =  &
    (/ 0.25000000000000000D0, 0.16666666666666667D0, &
       0.16666666666666667D0, 0.50000000000000000D0, &
       0.16666666666666667D0 /)

   REAL(KIND=dp), DIMENSION(5), PRIVATE :: STetra5P =  &
    (/-0.80000000000000000D0, 0.45000000000000000D0, &
       0.45000000000000000D0, 0.45000000000000000D0, &
       0.45000000000000000D0 /)
!------------------------------------------------------------------------------
! Tetrahedron - 11 point rule; exact integration of x^py^qz^r, p+q+r<=4
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UTetra11P =  &
    (/ 0.3247902050850455D+00, 0.4381969657060433D+00, &
       0.8992592373310454D-01, 0.1092714936292849D+00, &
       0.3389119319942253D-01, 0.5332363613904868D-01, &
       0.1935618747806815D+00, 0.4016250624424964D-01, &
       0.3878132182319405D+00, 0.7321489692875428D+00, &
       0.8066342495294049D-01 /)

   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VTetra11P =  &
    (/ 0.4573830181783998D+00, 0.9635325047480842D-01, &
       0.3499588148445295D+00, 0.1228957438582778D+00, &
       0.4736224692062527D-01, 0.4450376952468180D+00, &
       0.2165626476982170D+00, 0.8033385922433729D+00, &
       0.7030897281814283D-01, 0.1097836536360084D+00, &
       0.1018859284267242D-01 /)

   REAL(KIND=dp), DIMENSION(11), PRIVATE :: WTetra11P =  &
    (/ 0.1116787541193331D+00, 0.6966288385119494D-01, &
       0.5810783971325720D-01, 0.3424607753785182D+00, &
       0.7831772466208499D+00, 0.3688112094344830D+00, &
       0.5872345323698884D+00, 0.6178518963560731D-01, &
       0.4077342860913465D+00, 0.9607290317342082D-01, &
       0.8343823045787845D-01 /)

   REAL(KIND=dp), DIMENSION(11), PRIVATE :: STetra11P =  &
    (/ 0.1677896627448221D+00, 0.1128697325878004D+00, &
       0.1026246621329828D+00, 0.1583002576888426D+00, &
       0.3847841737508437D-01, 0.1061709382037234D+00, &
       0.5458124994014422D-01, 0.3684475128738168D-01, &
       0.1239234851349682D+00, 0.6832098141300300D-01, &
       0.3009586149124714D-01 /)
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! A SELECTION OF QUADRATURE RULES UP TO A HIGH ORDER DESIGNED FOR ECONOMIC
! INTEGRATION OF COMPLETE POLYNOMIALS
!
! These rules have been reproduced from the source "P. Solin, K. Segeth
! and I. Dolezel: Higher-Order Finite Element Methods",
! Chapman & Hall/CRC Press, 2003.
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Quadrilateral - 8-point rule for complete polynomials of order p<=5
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(8), PRIVATE :: UPQuad8 = &
       (/ 0.683130051063973d0,-0.683130051063973d0, 0.000000000000000d0, &
          0.000000000000000d0, 0.881917103688197d0, 0.881917103688197d0, &
          -0.881917103688197d0,-0.881917103688197d0 /)

   REAL(KIND=dp), DIMENSION(8), PRIVATE :: VPQuad8 = &
       (/ 0.000000000000000d0, 0.000000000000000d0, 0.683130051063973d0, &
         -0.683130051063973d0, 0.881917103688197d0,-0.881917103688197d0, &
          0.881917103688197d0,-0.881917103688197d0 /)

   REAL(KIND=dp), DIMENSION(8), PRIVATE :: SPQuad8 = &
       (/ 0.816326530612245d0, 0.816326530612245d0, 0.816326530612245d0, &
          0.816326530612245d0, 0.183673469387755d0, 0.183673469387755d0, &
          0.183673469387755d0, 0.183673469387755d0 /)
!------------------------------------------------------------------------------
! Quadrilateral - 12-point rule for complete polynomials of order p<=7
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(12), PRIVATE :: UPQuad12 = &
       (/ 0.925820099772551d0, -0.925820099772551d0,  0.000000000000000d0, &
          0.000000000000000d0,  0.805979782918599d0,  0.805979782918599d0, &
         -0.805979782918599d0, -0.805979782918599d0,  0.380554433208316d0, &
          0.380554433208316d0, -0.380554433208316d0, -0.380554433208316d0 /)

   REAL(KIND=dp), DIMENSION(12), PRIVATE :: VPQuad12 = &
       (/ 0.000000000000000d0,  0.000000000000000d0,  0.925820099772551d0, &
         -0.925820099772551d0,  0.805979782918599d0, -0.805979782918599d0, &
          0.805979782918599d0, -0.805979782918599d0,  0.380554433208316d0, &
         -0.380554433208316d0,  0.380554433208316d0, -0.380554433208316d0 /)

   REAL(KIND=dp), DIMENSION(12), PRIVATE :: SPQuad12 = &
       (/ 0.241975308641975d0, 0.241975308641975d0, 0.241975308641975d0, &
          0.241975308641975d0, 0.237431774690630d0, 0.237431774690630d0, &
          0.237431774690630d0, 0.237431774690630d0, 0.520592916667394d0, &
          0.520592916667394d0, 0.520592916667394d0, 0.520592916667394d0 /)
!------------------------------------------------------------------------------
! Quadrilateral - 20-point rule for complete polynomials of order p<=9
! Note: Some points lie outside the reference element
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(20), PRIVATE :: UPQuad20 = &
       (/ 0.112122576386656d1,  -0.112122576386656d1,   0.000000000000000d0,  &
          0.000000000000000d0,   0.451773049920657d0,  -0.451773049920657d0,  &
          0.000000000000000d0,   0.000000000000000d0,   0.891849420851512d0,  &
          0.891849420851512d0,  -0.891849420851512d0,  -0.891849420851512d0,  &
          0.824396370749276d0,   0.824396370749276d0,  -0.824396370749276d0,  &
         -0.824396370749276d0,   0.411623426336542d0,   0.411623426336542d0,  &
         -0.411623426336542d0,  -0.411623426336542d0 /)

   REAL(KIND=dp), DIMENSION(20), PRIVATE :: VPQuad20 = &
       (/ 0.000000000000000d0,  0.000000000000000d0,  0.112122576386656d1, &
         -0.112122576386656d1,  0.000000000000000d0,  0.000000000000000d0, &
          0.451773049920657d0, -0.451773049920657d0,  0.891849420851512d0, &
         -0.891849420851512d0,  0.891849420851512d0, -0.891849420851512d0, &
          0.411623426336542d0, -0.411623426336542d0,  0.411623426336542d0, &
         -0.411623426336542d0,  0.824396370749276d0, -0.824396370749276d0, &
         0.824396370749276d0, -0.824396370749276d0 /)

   REAL(KIND=dp), DIMENSION(20), PRIVATE :: SPQuad20 = &
       (/ 0.184758425074910d-1,  0.184758425074910d-1,  0.184758425074910d-1,  &
          0.184758425074910d-1,  0.390052939160735d0,   0.390052939160735d0,   &
          0.390052939160735d0,   0.390052939160735d0,   0.830951780264820d-1,  &
          0.830951780264820d-1,  0.830951780264820d-1,  0.830951780264820d-1,  &
          0.254188020152646d0,   0.254188020152646d0,   0.254188020152646d0,   &
          0.254188020152646d0,   0.254188020152646d0,   0.254188020152646d0,   &
          0.254188020152646d0,   0.254188020152646d0 /)
!------------------------------------------------------------------------------
! Quadrilateral - 25-point rule for complete polynomials of order p<=11
! Note: Some points lie outside the reference element
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(25), PRIVATE :: UPQuad25 = &
       (/0.000000000000000d0,   0.104440291540981d1,  -0.104440291540981d1,  &
         0.000000000000000d0,   0.000000000000000d0,   0.769799068396649d0,  &
        -0.769799068396649d0,   0.000000000000000d0,   0.000000000000000d0,  &
         0.935787012440540d0,   0.935787012440540d0,  -0.935787012440540d0,  &
        -0.935787012440540d0,   0.413491953449114d0,   0.413491953449114d0,  &
        -0.413491953449114d0,  -0.413491953449114d0,   0.883025508525690d0,  &
         0.883025508525690d0,  -0.883025508525690d0,  -0.883025508525690d0,  &
         0.575653595840465d0,   0.575653595840465d0,  -0.575653595840465d0,  &
        -0.575653595840465d0 /)

   REAL(KIND=dp), DIMENSION(25), PRIVATE :: VPQuad25 = &
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.000000000000000d0, &
          0.104440291540981d1,  -0.104440291540981d1,   0.000000000000000d0, &
          0.000000000000000d0,   0.769799068396649d0,  -0.769799068396649d0, &
          0.935787012440540d0,  -0.935787012440540d0,   0.935787012440540d0, &
         -0.935787012440540d0,   0.413491953449114d0,  -0.413491953449114d0, &
          0.413491953449114d0,  -0.413491953449114d0,   0.575653595840465d0, &
         -0.575653595840465d0,   0.575653595840465d0,  -0.575653595840465d0, &
          0.883025508525690d0,  -0.883025508525690d0,   0.883025508525690d0, &
          -0.883025508525690d0 /)

   REAL(KIND=dp), DIMENSION(25), PRIVATE :: SPQuad25 = &
       (/ 0.365379525585903d0,   0.277561655642040d-1,   0.277561655642040d-1, &
          0.277561655642040d-1,  0.277561655642040d-1,   0.244272057751754d0,  &
          0.244272057751754d0,   0.244272057751754d0,    0.244272057751754d0,  &
          0.342651038512290d-1,  0.342651038512290d-1,   0.342651038512290d-1, &
          0.342651038512290d-1,  0.308993036133713d0,    0.308993036133713d0,  &
          0.308993036133713d0,   0.308993036133713d0,    0.146684377651312d0,  &
          0.146684377651312d0,   0.146684377651312d0,    0.146684377651312d0,  &
          0.146684377651312d0,   0.146684377651312d0,    0.146684377651312d0,  &
          0.146684377651312d0  /)
!------------------------------------------------------------------------------
! Quadrilateral - 36-point rule for complete polynomials of order p<=13
! Note 1: Some points lie outside the reference element
! Note 2: Weights somewhat inaccurate?
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(36), PRIVATE :: UPQuad36 = &
       (/  0.108605615857397d1,  -0.108605615857397d1,   0.000000000000000d0, &
        0.000000000000000d0,   0.658208197042585d0,  -0.658208197042585d0,    &
        0.000000000000000d0,   0.000000000000000d0,   0.100130060299173d1,    &
        0.100130060299173d1,  -0.100130060299173d1,  -0.100130060299173d1,    &
        0.584636168775946d0,   0.584636168775946d0,  -0.584636168775946d0,    &
       -0.584636168775946d0,   0.246795612720261d0,   0.246795612720261d0,    &
       -0.246795612720261d0,  -0.246795612720261d0,   0.900258815287201d0,    &
        0.900258815287201d0,  -0.900258815287201d0,  -0.900258815287201d0,    &
        0.304720678579870d0,   0.304720678579870d0,  -0.304720678579870d0,    &
       -0.304720678579870d0,   0.929866705560780d0,   0.929866705560780d0,    &
       -0.929866705560780d0,  -0.929866705560780d0,   0.745052720131169d0,    &
        0.745052720131169d0,  -0.745052720131169d0,  -0.745052720131169d0  /)

   REAL(KIND=dp), DIMENSION(36), PRIVATE :: VPQuad36 = &
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.108605615857397d1,  &
       -0.108605615857397d1,   0.000000000000000d0,   0.000000000000000d0, &
        0.658208197042585d0,  -0.658208197042585d0,   0.100130060299173d1, &
       -0.100130060299173d1,   0.100130060299173d1,  -0.100130060299173d1, &
        0.584636168775946d0,  -0.584636168775946d0,   0.584636168775946d0, &
       -0.584636168775946d0,   0.246795612720261d0,  -0.246795612720261d0, &
        0.246795612720261d0,  -0.246795612720261d0,   0.304720678579870d0, &
       -0.304720678579870d0,   0.304720678579870d0,  -0.304720678579870d0, &
        0.900258815287201d0,  -0.900258815287201d0,   0.900258815287201d0, &
       -0.900258815287201d0,   0.745052720131169d0,  -0.745052720131169d0, &
        0.745052720131169d0,  -0.745052720131169d0,   0.929866705560780d0, &
       -0.929866705560780d0,   0.929866705560780d0,  -0.929866705560780d0 /)

   REAL(KIND=dp), DIMENSION(36), PRIVATE :: SPQuad36 = &
       (/ 0.565616969376400d-2,   0.565616969376400d-2,   0.565616969376400d-2,   &
       0.565616969376400d-2,  0.192443967470396d0,   0.192443967470396d0,  &
       0.192443967470396d0,   0.192443967470396d0,   0.516683297977300d-2, &
       0.516683297977300d-2,  0.516683297977300d-2,  0.516683297977300d-2, &
       0.200302559622138d0,   0.200302559622138d0,   0.200302559622138d0,  &
       0.200302559622138d0,   0.228125175912536d0,   0.228125175912536d0,  &
       0.228125175912536d0,   0.228125175912536d0,   0.117496926974491d0,  &
       0.117496926974491d0,   0.117496926974491d0,   0.117496926974491d0,  &
       0.117496926974491d0,   0.117496926974491d0,   0.117496926974491d0,  &
       0.117496926974491d0,   0.666557701862050d-1,  0.666557701862050d-1, &
       0.666557701862050d-1,  0.666557701862050d-1,  0.666557701862050d-1, &
       0.666557701862050d-1,  0.666557701862050d-1,  0.666557701862050d-1  /)
!------------------------------------------------------------------------------
! Quadrilateral - 45-point rule for complete polynomials of order p<=15
! Note: Some points lie outside the reference element
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(45), PRIVATE :: UPQuad45 = &
       (/ 0.000000000000000d0,   0.102731435771937d1,  -0.102731435771937d1,   &
        0.000000000000000d0,   0.000000000000000d0,   0.856766776147643d0,  &
       -0.856766776147643d0,   0.000000000000000d0,   0.000000000000000d0,  &
        0.327332998189723d0,  -0.327332998189723d0,   0.000000000000000d0,  &
        0.000000000000000d0,   0.967223740028505d0,   0.967223740028505d0,  &
       -0.967223740028505d0,  -0.967223740028505d0,   0.732168901749711d0,  &
        0.732168901749711d0,  -0.732168901749711d0,  -0.732168901749711d0,  &
        0.621974427996805d0,   0.621974427996805d0,  -0.621974427996805d0,  &
       -0.621974427996805d0,   0.321696694921009d0,   0.321696694921009d0,  &
       -0.321696694921009d0,  -0.321696694921009d0,   0.928618480068352d0,  &
        0.928618480068352d0,  -0.928618480068352d0,  -0.928618480068352d0,  &
        0.455124178121179d0,   0.455124178121179d0,  -0.455124178121179d0,  &
       -0.455124178121179d0,   0.960457474887516d0,   0.960457474887516d0,  &
       -0.960457474887516d0,  -0.960457474887516d0,   0.809863684081217d0,  &
        0.809863684081217d0,  -0.809863684081217d0,  -0.809863684081217d0   /)

   REAL(KIND=dp), DIMENSION(45), PRIVATE :: VPQuad45 = &
       (/ 0.000000000000000d0,   0.000000000000000d0,   0.000000000000000d0,   &
        0.102731435771937d1,  -0.102731435771937d1,   0.000000000000000d0, &
        0.000000000000000d0,   0.856766776147643d0,  -0.856766776147643d0, &
        0.000000000000000d0,   0.000000000000000d0,   0.327332998189723d0, &
       -0.327332998189723d0,   0.967223740028505d0,  -0.967223740028505d0, &
        0.967223740028505d0,  -0.967223740028505d0,   0.732168901749711d0, &
       -0.732168901749711d0,   0.732168901749711d0,  -0.732168901749711d0, &
        0.321696694921009d0,  -0.321696694921009d0,   0.321696694921009d0, &
       -0.321696694921009d0,   0.621974427996805d0,  -0.621974427996805d0, &
        0.621974427996805d0,  -0.621974427996805d0,   0.455124178121179d0, &
       -0.455124178121179d0,   0.455124178121179d0,  -0.455124178121179d0, &
        0.928618480068352d0,  -0.928618480068352d0,   0.928618480068352d0, &
       -0.928618480068352d0,   0.809863684081217d0,  -0.809863684081217d0, &
        0.809863684081217d0,  -0.809863684081217d0,   0.960457474887516d0, &
       -0.960457474887516d0,   0.960457474887516d0,  -0.960457474887516d0 /)

   REAL(KIND=dp), DIMENSION(45), PRIVATE :: SPQuad45 = &
       (/ -0.176897982720700d-2,   0.128167266175120d-1,   0.128167266175120d-1,   &
       0.128167266175120d-1,  0.128167266175120d-1,  0.119897873101347d0,   &
       0.119897873101347d0,   0.119897873101347d0,   0.119897873101347d0,   &
       0.210885452208801d0,   0.210885452208801d0,   0.210885452208801d0,   &
       0.210885452208801d0,   0.639272012821500d-2,  0.639272012821500d-2,  &
       0.639272012821500d-2,  0.639272012821500d-2,  0.104415680788580d0,   &
       0.104415680788580d0,   0.104415680788580d0,   0.104415680788580d0,   &
       0.168053047203816d0,   0.168053047203816d0,   0.168053047203816d0,   &
       0.168053047203816d0,   0.168053047203816d0,   0.168053047203816d0,   &
       0.168053047203816d0,   0.168053047203816d0,   0.761696944522940d-1,  &
       0.761696944522940d-1,  0.761696944522940d-1,  0.761696944522940d-1,  &
       0.761696944522940d-1,  0.761696944522940d-1,  0.761696944522940d-1,  &
       0.761696944522940d-1,  0.287941544000640d-1,  0.287941544000640d-1,  &
       0.287941544000640d-1,  0.287941544000640d-1,  0.287941544000640d-1,  &
       0.287941544000640d-1,  0.287941544000640d-1,  0.287941544000640d-1 /)
!------------------------------------------------------------------------------
! Quadrilateral - 60-point rule for complete polynomials of order p<=17
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(60), PRIVATE :: UPQuad60 = &
       (/ 0.989353074512600d0,  -0.989353074512600d0,   0.000000000000000d0,   &
        0.000000000000000d0,   0.376285207157973d0,  -0.376285207157973d0,  &
        0.000000000000000d0,   0.000000000000000d0,   0.978848279262233d0,  &
        0.978848279262233d0,  -0.978848279262233d0,  -0.978848279262233d0,  &
        0.885794729164116d0,   0.885794729164116d0,  -0.885794729164116d0,  &
       -0.885794729164116d0,   0.171756123838348d0,   0.171756123838348d0,  &
       -0.171756123838348d0,  -0.171756123838348d0,   0.590499273806002d0,  &
        0.590499273806002d0,  -0.590499273806002d0,  -0.590499273806002d0,  &
        0.319505036634574d0,   0.319505036634574d0,  -0.319505036634574d0,  &
       -0.319505036634574d0,   0.799079131916863d0,   0.799079131916863d0,  &
       -0.799079131916863d0,  -0.799079131916863d0,   0.597972451929457d0,  &
        0.597972451929457d0,  -0.597972451929457d0,  -0.597972451929457d0,  &
        0.803743962958745d0,   0.803743962958745d0,  -0.803743962958745d0,  &
       -0.803743962958745d0,   0.583444817765510d-1,  0.583444817765510d-1, &
       -0.583444817765510d-1, -0.583444817765510d-1,  0.936506276127495d0,  &
        0.936506276127495d0,  -0.936506276127495d0,  -0.936506276127495d0,  &
        0.347386316166203d0,   0.347386316166203d0,  -0.347386316166203d0,  &
       -0.347386316166203d0,   0.981321179805452d0,   0.981321179805452d0,  &
       -0.981321179805452d0,  -0.981321179805452d0,   0.706000287798646d0,  &
        0.706000287798646d0,  -0.706000287798646d0,  -0.706000287798646d0 /)

   REAL(KIND=dp), DIMENSION(60), PRIVATE :: VPQuad60 = &
       (/0.000000000000000d0,   0.000000000000000d0,   0.989353074512600d0,  &
       -0.989353074512600d0,   0.000000000000000d0,   0.000000000000000d0,  &
        0.376285207157973d0,  -0.376285207157973d0,   0.978848279262233d0,  &
       -0.978848279262233d0,   0.978848279262233d0,  -0.978848279262233d0,  &
        0.885794729164116d0,  -0.885794729164116d0,   0.885794729164116d0,  &
       -0.885794729164116d0,   0.171756123838348d0,  -0.171756123838348d0,  &
        0.171756123838348d0,  -0.171756123838348d0,   0.319505036634574d0,  &
       -0.319505036634574d0,   0.319505036634574d0,  -0.319505036634574d0,  &
        0.590499273806002d0,  -0.590499273806002d0,   0.590499273806002d0,  &
       -0.590499273806002d0,   0.597972451929457d0,  -0.597972451929457d0,  &
        0.597972451929457d0,  -0.597972451929457d0,   0.799079131916863d0,  &
       -0.799079131916863d0,   0.799079131916863d0,  -0.799079131916863d0,  &
        0.583444817765510d-1, -0.583444817765510d-1,  0.583444817765510d-1, &
       -0.583444817765510d-1,  0.803743962958745d0,  -0.803743962958745d0,  &
        0.803743962958745d0,  -0.803743962958745d0,   0.347386316166203d0,  &
       -0.347386316166203d0,   0.347386316166203d0,  -0.347386316166203d0,  &
        0.936506276127495d0,  -0.936506276127495d0,   0.936506276127495d0,  &
       -0.936506276127495d0,   0.706000287798646d0,  -0.706000287798646d0,  &
        0.706000287798646d0,  -0.706000287798646d0,   0.981321179805452d0,  &
       -0.981321179805452d0,   0.981321179805452d0,  -0.981321179805452d0  /)

   REAL(KIND=dp), DIMENSION(60), PRIVATE :: SPQuad60 = &
       (/ 0.206149159199910d-1,   0.206149159199910d-1,   0.206149159199910d-1,   &
       0.206149159199910d-1,  0.128025716179910d0,    0.128025716179910d0,  &
       0.128025716179910d0,   0.128025716179910d0,    0.551173953403200d-2, &
       0.551173953403200d-2,  0.551173953403200d-2,   0.551173953403200d-2, &
       0.392077124571420d-1,  0.392077124571420d-1,   0.392077124571420d-1, &
       0.392077124571420d-1,  0.763969450798630d-1,   0.763969450798630d-1, &
       0.763969450798630d-1,  0.763969450798630d-1,   0.141513729949972d0,  &
       0.141513729949972d0,   0.141513729949972d0,    0.141513729949972d0,  &
       0.141513729949972d0,   0.141513729949972d0,    0.141513729949972d0,  &
       0.141513729949972d0,   0.839032793637980d-1,   0.839032793637980d-1, &
       0.839032793637980d-1,  0.839032793637980d-1,   0.839032793637980d-1, &
       0.839032793637980d-1,  0.839032793637980d-1,   0.839032793637980d-1, &
       0.603941636496850d-1,  0.603941636496850d-1,   0.603941636496850d-1, &
       0.603941636496850d-1,  0.603941636496850d-1,   0.603941636496850d-1, &
       0.603941636496850d-1,  0.603941636496850d-1,   0.573877529692130d-1, &
       0.573877529692130d-1,  0.573877529692130d-1,   0.573877529692130d-1, &
       0.573877529692130d-1,  0.573877529692130d-1,   0.573877529692130d-1, &
       0.573877529692130d-1,  0.219225594818640d-1,   0.219225594818640d-1, &
       0.219225594818640d-1,  0.219225594818640d-1,   0.219225594818640d-1, &
       0.219225594818640d-1,  0.219225594818640d-1,   0.219225594818640d-1  /)


!------------------------------------------------------------------------------
! These rules have been reproduced from the paper
! Ethan J. Kubatko, Benjamin A. Yeager, Ashley L. Magg and I. Dolezel:
! "New computationally efficient quadrature formulas for triangular
! prism elements.", Computers & Fluids 73 (2013) 187–201.
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
! Wedge - 4 point rule, 2nd order polynoms
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: SWedge4P = [ &
       1.000000000000000d0, &
       1.000000000000000d0, &
       1.000000000000000d0, &
       1.000000000000000d0]
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: UWedge4P = [ &
       0.483163247594393d0, &
       -0.605498860309242d0, &
       -0.605498860309242d0, &
       -0.605498860309242d0]
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: VWedge4P = [ &
       -0.741581623797196d0, &
       0.469416096821288d0, &
       -0.530583903178712d0, &
       -0.530583903178712d0]
   REAL(KIND=dp), DIMENSION(4), PRIVATE :: WWedge4P = [ &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.816496580927726d0, &
       -0.816496580927726d0 ]

!------------------------------------------------------------------------------
! Wedge - n=5, p=2
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: SWedge5P = [ &
       0.333333333333333d0, &
       0.333333333333333d0, &
       0.333333333333333d0, &
       1.500000000000000d0, &
       1.500000000000000d0 ]
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: UWedge5P = [ &
       -1.000000000000000d0, &
       -1.000000000000000d0, &
       1.000000000000000d0, &
       -0.333333333333333d0, &
       -0.333333333333333d0 ]
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: VWedge5P = [ &
       -1.000000000000000d0, &
       1.000000000000000d0, &
       -1.000000000000000d0, &
       -0.333333333333333d0, &
       -0.333333333333333d0 ]
   REAL(KIND=dp), DIMENSION(5), PRIVATE :: WWedge5P = [ &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.666666666666667d0, &
       -0.666666666666667d0 ]

!------------------------------------------------------------------------------
! Wedge - n=6, p=3
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: SWedge6P = [ &
       0.742534890852309d0, &
       0.375143463443327d0, &
       0.495419047908462d0, &
       0.523999970843238d0, &
       0.980905839025611d0, &
       0.881996787927053d0 ]
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: UWedge6P = [ &
       0.240692796349625d0, &
       -0.968326281451138d0, &
       0.467917833640195d0, &
       -0.786144119530819d0, &
       -0.484844897886675d0, &
       -0.559053711932125d0 ]
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: VWedge6P = [ &
       -0.771991660873412d0, &
       -0.568046512457875d0, &
       -0.549342790168347d0, &
       0.362655041695561d0, &
       -0.707931130717342d0, &
       0.260243325246813d0 ]
   REAL(KIND=dp), DIMENSION(6), PRIVATE :: WWedge6P = [ &
       0.614747128207527d0, &
       0.676689529541421d0, &
       -0.599905857322635d0, &
       -0.677609795694786d0, &
       -0.502482717716373d0, &
       0.493010512161538d0 ]

!------------------------------------------------------------------------------
! Wedge - n=7, p=3
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: SWedge7P = [ &
       -2.250000000000000d0, &
       1.041666666666667d0, &
       1.041666666666667d0, &
       1.041666666666667d0, &
       1.041666666666667d0, &
       1.041666666666667d0, &
       1.041666666666667d0 ]
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: UWedge7P = [ &
       -0.333333333333333d0, &
       -0.600000000000000d0, &
       -0.600000000000000d0, &
       0.200000000000000d0, &
       -0.600000000000000d0, &
       -0.600000000000000d0, &
       0.200000000000000d0 ]
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: VWedge7P = [ &
       -0.333333333333333d0, &
       -0.600000000000000d0, &
       0.200000000000000d0, &
       -0.600000000000000d0, &
       -0.600000000000000d0, &
       0.200000000000000d0, &
       -0.600000000000000d0 ]
   REAL(KIND=dp), DIMENSION(7), PRIVATE :: WWedge7P = [ &
       0.000000000000000d0, &
       0.461880215351701d0, &
       0.461880215351701d0, &
       0.461880215351701d0, &
       -0.461880215351701d0, &
       -0.461880215351701d0, &
       -0.461880215351701d0 ]

!------------------------------------------------------------------------------
! Wedge - n=10, p=4
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: SWedge10P = [ &
       0.111155943811228d0, &
       0.309060899887509d0, &
       0.516646862442958d0, &
       0.567975205132714d0, &
       0.382742555939017d0, &
       0.355960928492268d0, &
       0.108183228294342d0, &
       0.126355242780924d0, &
       0.587370828592853d0, &
       0.934548304626188d0 ]
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: UWedge10P = [ &
       0.812075900047562d0, &
       -0.792166223585545d0, &
       -0.756726179789306d0, &
       -0.552495167978340d0, &
       -0.357230019521233d0, &
       -0.987225392999058d0, &
       -0.816603728785918d0, &
       0.423489172633859d0, &
       0.363041084609230d0, &
       -0.175780343149613d0 ]
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: VWedge10P = [ &
       -0.986242751499303d0, &
       0.687201105597868d0, &
       -0.731311840596107d0, &
       0.015073398439985d0, &
       0.126888850505978d0, &
       0.082647545710800d0, &
       -0.915066171481315d0, &
       -1.112801167237130d0, &
       -0.499011410082669d0, &
       -0.654971142379686d0 ]
   REAL(KIND=dp), DIMENSION(10), PRIVATE :: WWedge10P = [ &
       0.850716248413834d0, &
       -0.115214772515700d0, &
       -0.451491675441927d0, &
       -0.824457000064439d0, &
       0.855349689995606d0, &
       0.452976444667786d0, &
       0.997939285245240d0, &
       -0.963298774205756d0, &
       -0.299892769705443d0, &
       0.367947041936472d0 ]

!------------------------------------------------------------------------------
! Wedge - n=11, p=4
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: SWedge11P = [ &
       0.545658450421913d0, &
       0.545658450421913d0, &
       0.545658450421913d0, &
       0.431647899262139d0, &
       0.249954808368331d0, &
       0.249954808368331d0, &
       0.249954808368331d0, &
       0.431647899262139d0, &
       0.249954808368331d0, &
       0.249954808368331d0, &
       0.249954808368331d0 ]
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: UWedge11P = [ &
       -0.062688380276010d0, &
       -0.062688380276010d0, &
       -0.874623239447980d0, &
       -0.333333333333333d0, &
       -0.798519188402179d0, &
       -0.798519188402179d0, &
       0.597038376804357d0, &
       -0.333333333333333d0, &
       -0.798519188402179d0, &
       -0.798519188402179d0, &
       0.597038376804357d0 ]
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: VWedge11P = [ &
       -0.062688380276010d0, &
       -0.874623239447980d0, &
       -0.062688380276010d0, &
       -0.333333333333333d0, &
       -0.798519188402179d0, &
       0.597038376804357d0, &
       -0.798519188402179d0, &
       -0.333333333333333d0, &
       -0.798519188402179d0, &
       0.597038376804357d0, &
       -0.798519188402179d0 ]
   REAL(KIND=dp), DIMENSION(11), PRIVATE :: WWedge11P = [ &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.866861974009030d0, &
       0.675639823682265d0, &
       0.675639823682265d0, &
       0.675639823682265d0, &
       -0.866861974009030d0, &
       -0.675639823682265d0, &
       -0.675639823682265d0, &
       -0.675639823682265d0 ]

!------------------------------------------------------------------------------
! Wedge - n=14, p=5
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: SWedge14P = [ &
       0.087576186678730d0, &
       0.229447629454892d0, &
       0.229447629454891d0, &
       0.833056798542985d0, &
       0.166443428304729d0, &
       0.166443428304729d0, &
       0.376993270712316d0, &
       0.170410864470884d0, &
       0.298194157223163d0, &
       0.298194157223162d0, &
       0.376993270712316d0, &
       0.170410864470884d0, &
       0.298194157223163d0, &
       0.298194157223162d0 ]
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: UWedge14P = [ &
       -0.955901336867574d0, &
       -0.051621305926029d0, &
       -1.017063354038640d0, &
       -0.297388746460523d0, &
       0.774849157169622d0, &
       -0.849021640062097d0, &
       -0.665292008657551d0, &
       -0.012171148087346d0, &
       -0.734122164680096d0, &
       0.334122164680066d0, &
       -0.665292008657551d0, &
       -0.012171148087346d0, &
       -0.734122164680096d0, &
       0.334122164680066d0 ]
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: VWedge14P = [ &
       -0.955901336867577d0, &
       -1.017063354038640d0, &
       -0.051621305926032d0, &
       -0.297388746460521d0, &
       -0.849021640062096d0, &
       0.774849157169620d0, &
       -0.665292008657551d0, &
       -0.012171148087349d0, &
       0.334122164680066d0, &
       -0.734122164680096d0, &
       -0.665292008657551d0, &
       -0.012171148087349d0, &
       0.334122164680066d0, &
       -0.734122164680096d0 ]
   REAL(KIND=dp), DIMENSION(14), PRIVATE :: WWedge14P = [ &
       0.000000000000286d0, &
       0.000000000000038d0, &
       0.000000000000038d0, &
       0.000000000000008d0, &
       0.000000000000080d0, &
       0.000000000000080d0, &
       0.756615409654429d0, &
       0.600149379161583d0, &
       0.808115770496521d0, &
       0.808115770496521d0, &
       -0.756615409654429d0, &
       -0.600149379161583d0, &
       -0.808115770496521d0, &
       -0.808115770496521d0 ]

!------------------------------------------------------------------------------
! Wedge - n=15, p=5
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: SWedge15P = [ &
       0.213895020288765d0, &
       0.141917375616806d0, &
       0.295568859378071d0, &
       0.256991945593379d0, &
       0.122121979248381d0, &
       0.175194917962627d0, &
       0.284969106392719d0, &
       0.323336131783334d0, &
       0.159056110329063d0, &
       0.748067388709644d0, &
       0.280551223607231d0, &
       0.147734016552639d0, &
       0.259874920383688d0, &
       0.235144061421191d0, &
       0.355576942732463d0 ]
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: UWedge15P = [ &
       -0.820754415297359d0, &
       0.611831616907812d0, &
       -0.951379065092975d0, &
       0.200535109198601d0, &
       -0.909622841605196d0, &
       0.411514133287729d0, &
       -0.127534496411879d0, &
       -0.555217727817199d0, &
       0.706942532529193d0, &
       -0.278092963133809d0, &
       -0.057824844208300d0, &
       -0.043308436222116d0, &
       -0.774478920734726d0, &
       -0.765638443571458d0, &
       -0.732830649614460d0 ]
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: VWedge15P = [ &
       0.701020947925133d0, &
       -0.869995576950389d0, &
       -0.087091980055873d0, &
       -0.783721735474016d0, &
       -0.890218158063352d0, &
       -0.725374126531787d0, &
       -0.953467283619037d0, &
       -0.530472194607007d0, &
       -0.782248553944540d0, &
       -0.291936530517119d0, &
       -0.056757587543798d0, &
       0.012890722780611d0, &
       0.476188541042454d0, &
       0.177195164202219d0, &
       -0.737447982744191d0 ]
   REAL(KIND=dp), DIMENSION(15), PRIVATE :: WWedge15P = [ &
       -0.300763696502910d0, &
       0.348546607420888d0, &
       0.150656042323906d0, &
       -0.844285153176719d0, &
       0.477120081549168d0, &
       0.864653509536562d0, &
       0.216019762875977d0, &
       0.873409672725819d0, &
       -0.390653804976705d0, &
       -0.126030507204870d0, &
       0.539907869785112d0, &
       -0.776314479909204d0, &
       0.745875967497062d0, &
       -0.888355356215127d0, &
       -0.651653242952189d0 ]

!------------------------------------------------------------------------------
! Wedge - n=16, p=5
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: SWedge16P = [ &
       0.711455555931488d0, &
       0.224710067228267d0, &
       0.224710067228267d0, &
       0.224710067228267d0, &
       0.185661421316158d0, &
       0.185661421316158d0, &
       0.185661421316158d0, &
       0.250074285747794d0, &
       0.250074285747794d0, &
       0.250074285747794d0, &
       0.185661421316158d0, &
       0.185661421316158d0, &
       0.185661421316158d0, &
       0.250074285747794d0, &
       0.250074285747794d0, &
       0.250074285747794d0 ]
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: UWedge16P = [ &
       -0.333333333333333d0, &
       -0.025400070899509d0, &
       -0.025400070899509d0, &
       -0.949199858200983d0, &
       -0.108803790659256d0, &
       -0.108803790659256d0, &
       -0.782392418681488d0, &
       -0.798282108034583d0, &
       -0.798282108034583d0, &
       0.596564216069166d0, &
       -0.108803790659256d0, &
       -0.108803790659256d0, &
       -0.782392418681488d0, &
       -0.798282108034583d0, &
       -0.798282108034583d0, &
       0.596564216069166d0 ]
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: VWedge16P = [ &
       -0.333333333333333d0, &
       -0.025400070899509d0, &
       -0.949199858200983d0, &
       -0.025400070899509d0, &
       -0.108803790659256d0, &
       -0.782392418681488d0, &
       -0.108803790659256d0, &
       -0.798282108034583d0, &
       0.596564216069166d0, &
       -0.798282108034583d0, &
       -0.108803790659256d0, &
       -0.782392418681488d0, &
       -0.108803790659256d0, &
       -0.798282108034583d0, &
       0.596564216069166d0, &
       -0.798282108034583d0 ]
   REAL(KIND=dp), DIMENSION(16), PRIVATE :: WWedge16P = [ &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.871002934865444d0, &
       0.871002934865444d0, &
       0.871002934865444d0, &
       0.570426980705159d0, &
       0.570426980705159d0, &
       0.570426980705159d0, &
       -0.871002934865444d0, &
       -0.871002934865444d0, &
       -0.871002934865444d0, &
       -0.570426980705159d0, &
       -0.570426980705159d0, &
       -0.570426980705159d0 ]

!------------------------------------------------------------------------------
! Wedge - n=24, p=6
!------------------------------------------------------------------------------
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: SWedge24P = [ &
       0.168480079940386d0, &
       0.168480079940386d0, &
       0.168480079940386d0, &
       0.168480079940386d0, &
       0.168480079940386d0, &
       0.168480079940386d0, &
       0.000079282874851d0, &
       0.000079282874851d0, &
       0.000079282874851d0, &
       0.544286440652304d0, &
       0.026293733850586d0, &
       0.283472344926041d0, &
       0.283472344926041d0, &
       0.283472344926041d0, &
       0.115195615637235d0, &
       0.115195615637235d0, &
       0.115195615637235d0, &
       0.026293733850586d0, &
       0.283472344926041d0, &
       0.283472344926041d0, &
       0.283472344926041d0, &
       0.115195615637235d0, &
       0.115195615637235d0, &
       0.115195615637235d0 ]
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: UWedge24P = [ &
       0.513019949700545d0, &
       -0.582925557762337d0, &
       -0.582925557762337d0, &
       0.513019949700545d0, &
       -0.930094391938207d0, &
       -0.930094391938207d0, &
       -1.830988812620400d0, &
       -1.830988812620400d0, &
       2.661977625240810d0, &
       -0.333333333333333d0, &
       -0.333333333333333d0, &
       -0.098283514203544d0, &
       -0.098283514203544d0, &
       -0.803432971592912d0, &
       -0.812603471654584d0, &
       -0.812603471654584d0, &
       0.625206943309168d0, &
       -0.333333333333333d0, &
       -0.098283514203544d0, &
       -0.098283514203544d0, &
       -0.803432971592912d0, &
       -0.812603471654584d0, &
       -0.812603471654584d0, &
       0.625206943309168d0 ]
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: VWedge24P = [ &
       -0.930094391938207d0, &
       -0.930094391938207d0, &
       0.513019949700545d0, &
       -0.582925557762337d0, &
       -0.582925557762337d0, &
       0.513019949700545d0, &
       -1.830988812620400d0, &
       2.661977625240810d0, &
       -1.830988812620400d0, &
       -0.333333333333333d0, &
       -0.333333333333333d0, &
       -0.098283514203544d0, &
       -0.803432971592912d0, &
       -0.098283514203544d0, &
       -0.812603471654584d0, &
       0.625206943309168d0, &
       -0.812603471654584d0, &
       -0.333333333333333d0, &
       -0.098283514203544d0, &
       -0.803432971592912d0, &
       -0.098283514203544d0, &
       -0.812603471654584d0, &
       0.625206943309168d0, &
       -0.812603471654584d0 ]
   REAL(KIND=dp), DIMENSION(24), PRIVATE :: WWedge24P = [ &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       0.000000000000000d0, &
       1.250521622121900d0, &
       0.685008566774710d0, &
       0.685008566774710d0, &
       0.685008566774710d0, &
       0.809574716992997d0, &
       0.809574716992997d0, &
       0.809574716992997d0, &
       -1.250521622121900d0, &
       -0.685008566774710d0, &
       -0.685008566774710d0, &
       -0.685008566774710d0, &
       -0.809574716992997d0, &
       -0.809574716992997d0, &
       -0.809574716992997d0 ]

! continue wedge rules here
!------------------------------------------------------------------------------
! Wedge - n=27, p=6
!------------------------------------------------------------------------------
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: SWedge27P = [ &
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: UWedge27P = [ &
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: VWedge27P = [ &
!   REAL(KIND=dp), DIMENSION(27), PRIVATE :: WWedge27P = [ &



 CONTAINS


!------------------------------------------------------------------------------
!> Subroutine to compute Fejer integration points by FFT.
!------------------------------------------------------------------------------
   SUBROUTINE ComputeFejerPoints1D( Points,Weights,n )
!------------------------------------------------------------------------------
    INTEGER :: n
    REAL(KIND=dp) :: Points(n),Weights(n)
!------------------------------------------------------------------------------
    INTEGER :: i,j,k,l,m
    TYPE(C_FFTCMPLX) :: W(n+1)
    REAL(KIND=dp) :: arr((n+1)/2+1), V(n+2)

    DO i=1,(n+1)/2
      Points(i) = COS(i*PI/(n+1._dp))
      Points(n-i+1) = -Points(i)
    END DO

    k = 0
    DO i=1,n+1,2
      k = k + 1
      arr(k) = i
    END DO

    V = 0._dp
    DO i=1,k
      V(i) = 2._dp/(arr(i)*(arr(i)-2))
    END DO
    V(k+1) = 1._dp/arr(k)

    DO i=1,n+1
      W(i) % rval = -(V(i)+V(n-i+3))
    END DO
    CALL frfftb(n+1,W,V)

    Weights(1:n) = V(2:n+1)/(n+1)
    DO i=1,n/2
      Weights(i) = (Weights(i) + Weights(n-i+1))/2
      Weights(n-i+1) = Weights(i)
    END DO
    Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n))
!------------------------------------------------------------------------------
  END SUBROUTINE ComputeFejerPoints1D
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> Subroutine to compute Gaussian integration points and weights in [-1,1]
!> as roots of Legendre polynomials.
!------------------------------------------------------------------------------
   SUBROUTINE ComputeGaussPoints1D( Points,Weights,n )
!------------------------------------------------------------------------------
     INTEGER :: n
     REAL(KIND=dp) :: Points(n),Weights(n)
!------------------------------------------------------------------------------
     REAL(KIND=dp)   :: A(n/2,n/2),s,x,Work(8*n)
     COMPLEX(KIND=dp) :: Eigs(n/2)
     REAL(KIND=dp)   :: P(n+1),Q(n),P0(n),P1(n+1)
     INTEGER :: i,j,k,np,info
!------------------------------------------------------------------------------
! One point is trivial
!------------------------------------------------------------------------------
     IF ( n <= 1 ) THEN
       Points(1)  = 0.0d0
       Weights(1) = 2.0d0
       RETURN
     END IF
!------------------------------------------------------------------------------
! Compute coefficients of n:th Legendre polynomial from the recurrence:
!
! (i+1)P_{i+1}(x) = (2i+1)*x*P_i(x) - i*P_{i-1}(x), P_{0} = 1; P_{1} = x;
!
! CAVEAT: Computed coefficients inaccurate for n > ~15
!------------------------------------------------------------------------------
     P = 0
     P0(1) = 1
     P1(1) = 1
     P1(2) = 0

     DO i=1,n-1
       P(1:i+1) = (2*i+1) * P1(1:i+1)  / (i+1)
       P(3:i+2) = P(3:i+2) - i*P0(1:i) / (i+1)
       P0(1:i+1) = P1(1:i+1); P1(1:i+2) = P(1:i+2)
     END DO
!------------------------------------------------------------------------------
! Odd n implicates zero as one of the roots...
!------------------------------------------------------------------------------
     np = n - MOD(n,2)
!------------------------------------------------------------------------------
!  Variable substitution: y=x^2
!------------------------------------------------------------------------------
     np = np / 2
     DO i=1,np+1
       P(i) = P(2*i-1)
     END DO
!------------------------------------------------------------------------------
! Solve the roots of the polynomial by forming a matrix whose characteristic
! polynomial is the n:th Legendre polynomial and solving for the eigenvalues.
! Dunno if this is a very good method....
!------------------------------------------------------------------------------
     A=0
     DO i=1,np-1
       A(i,i+1) = 1
     END DO

     DO i=1,np
       A(np,i) = -P(np+2-i) / P(1)
     END DO
     CALL DGEEV( 'N','N',np,A,n/2,Points,P0,Work,1,Work,1,Work,8*n,info )
!    CALL EigenValues( A, np, Eigs )
!    Points(1:np) = REAL( Eigs(1:np) )

!------------------------------------------------------------------------------
! Backsubstitute from y=x^2
!------------------------------------------------------------------------------
     Q(1:np+1) = P(1:np+1)
     P = 0
     DO i=1,np+1
       P(2*i-1) = Q(i)
     END DO

     Q(1:np) = Points(1:np)
     DO i=1,np
       Points(2*i-1) = +SQRT( Q(i) )
       Points(2*i)   = -SQRT( Q(i) )
     END DO
     IF ( MOD(n,2) == 1 ) Points(n) = 0.0d0

     CALL DerivPoly( n,Q,P )
     CALL RefineRoots( n,P,Q,Points )
!------------------------------------------------------------------------------
! Check for roots
!------------------------------------------------------------------------------
     DO i=1,n
       s = EvalPoly( n,P,Points(i) )
       IF ( ABS(s) > 1.0d-12 ) THEN
         CALL Warn( 'ComputeGaussPoints1D', &
                 '-------------------------------------------------------' )
         CALL Warn( 'ComputeGaussPoints1D', 'Computed integration point' )
         CALL Warn( 'ComputeGaussPoints1D', 'seems to be inaccurate: ' )
         WRITE( Message, * ) 'Points req.: ',n
         CALL Warn( 'ComputeGaussPoints1D', Message )
         WRITE( Message, * ) 'Residual: ',s
         CALL Warn( 'ComputeGaussPoints1D', Message )
         WRITE( Message, * ) 'Point: +-', SQRT(Points(i))
         CALL Warn( 'ComputeGaussPoints1D', Message )
         CALL Warn( 'ComputeGaussPoints1D', &
                 '-------------------------------------------------------' )
       END IF
     END DO
!------------------------------------------------------------------------------
! Finally, the integration weights equal to
!
! W_i = 2/( (1-x_i^2)*Q(x_i)^2 ), x_i is the i:th root, and Q(x) = dP(x) / dx
!------------------------------------------------------------------------------
     CALL DerivPoly( n,Q,P )

     DO i=1,n
       s = EvalPoly( n-1,Q,Points(i) )
       Weights(i) = 2 / ((1-Points(i)**2)*s**2);
     END DO
!------------------------------------------------------------------------------
! ... make really sure the weights add up:
!------------------------------------------------------------------------------
     Weights(1:n) = 2 * Weights(1:n) / SUM(Weights(1:n))

CONTAINS

!--------------------------------------------------------------------------

   FUNCTION EvalPoly( n,P,x ) RESULT(s)
     INTEGER :: i,n
     REAL(KIND=dp) :: P(n+1),x,s

     s = 0.0d0
     DO i=1,n+1
       s = s * x + P(i)
     END DO
   END FUNCTION EvalPoly

!--------------------------------------------------------------------------

   SUBROUTINE DerivPoly( n,Q,P )
     INTEGER :: i,n
     REAL(KIND=dp) :: Q(n),P(n+1)

     DO i=1,n
       Q(i) = P(i)*(n-i+1)
     END DO
   END SUBROUTINE DerivPoly

!--------------------------------------------------------------------------

   SUBROUTINE RefineRoots( n,P,Q,Points )
     INTEGER :: i,j,n
     REAL(KIND=dp) :: P(n+1),Q(n),Points(n)

     REAL(KIND=dp) :: x,s
     INTEGER, PARAMETER :: MaxIter = 100

     DO i=1,n
       x = Points(i)
       DO j=1,MaxIter
         s = EvalPoly(n,P,x) / EvalPoly(n-1,Q,x)
         x = x - s
         IF ( ABS(s) <= ABS(x)*EPSILON(s) ) EXIT
       END DO
       IF ( ABS(EvalPoly(n,P,x))<ABS(EvalPoly(n,P,Points(i))) ) THEN
         IF ( ABS(x-Points(i))<1.0d-8 ) Points(i) = x
       END IF
     END DO
   END SUBROUTINE RefineRoots

!--------------------------------------------------------------------------
 END SUBROUTINE ComputeGaussPoints1D
!--------------------------------------------------------------------------

!--------------------------------------------------------------------------
   FUNCTION GaussPointsInitialized() RESULT(g)
!--------------------------------------------------------------------------
     LOGICAL :: g
!--------------------------------------------------------------------------
     g=Ginit
!--------------------------------------------------------------------------
   END FUNCTION GaussPointsInitialized
!--------------------------------------------------------------------------

! !------------------------------------------------------------------------------
!    SUBROUTINE GaussPointsInit
! !------------------------------------------------------------------------------
!      INTEGER :: i,n,istat
!      INTEGER :: nthreads, thread, omp_get_max_threads, omp_get_thread_num
!
!      IF ( .NOT. ALLOCATED(IntegStuff) ) THEN
! !$omp barrier
! !$omp critical
!        IF ( .NOT. GInit ) THEN
!          GInit = .TRUE.
!          nthreads = 1
! !$       nthreads = omp_get_max_threads()
!          ALLOCATE( IntegStuff(nthreads) )
!          DO i=1,nthreads
!            IntegStuff(i) % u => NULL()
!            IntegStuff(i) % v => NULL()
!            IntegStuff(i) % w => NULL()
!            IntegStuff(i) % s => NULL()
!          END DO
!          DO n=1,MAXN
!            CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n )
!          END DO
!        END IF
! !$omp end critical
!      END IF
!
!      thread = 1
! !$   thread = omp_get_thread_num()+1
!     ALLOCATE( IntegStuff(thread) % u(MAX_INTEGRATION_POINTS), &
!               IntegStuff(thread) % v(MAX_INTEGRATION_POINTS), &
!               IntegStuff(thread) % w(MAX_INTEGRATION_POINTS), &
!               IntegStuff(thread) % s(MAX_INTEGRATION_POINTS), STAT=istat )
!
!     IF ( istat /= 0 ) THEN
!       CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' )
!       STOP
!     END IF
! !------------------------------------------------------------------------------
! !  END SUBROUTINE GaussPointsInit
! !------------------------------------------------------------------------------

!------------------------------------------------------------------------------
   SUBROUTINE GaussPointsInit
!------------------------------------------------------------------------------
     INTEGER :: i,n,istat

     IF ( .NOT. GInit ) THEN
        DO n=1,MAXN
          CALL ComputeGaussPoints1D( Points(1:n,n),Weights(1:n,n),n )
        END DO
        GInit = .TRUE.
     END IF

     ALLOCATE( IntegStuff % u(MAX_INTEGRATION_POINTS), &
               IntegStuff % v(MAX_INTEGRATION_POINTS), &
               IntegStuff % w(MAX_INTEGRATION_POINTS), &
               IntegStuff % s(MAX_INTEGRATION_POINTS), STAT=istat )
     IntegStuff % u = 0._dp
     IntegStuff % v = 0._dp
     IntegStuff % w = 0._dp
     IntegStuff % s = 0._dp
     IF ( istat /= 0 ) THEN
       CALL Fatal( 'GaussPointsInit', 'Memory allocation error.' )
     END IF
!------------------------------------------------------------------------------
  END SUBROUTINE GaussPointsInit
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
   FUNCTION GaussPoints0D( n ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: n
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!     INTEGER :: thread, omp_get_thread_num

      IF ( .NOT. GInit ) CALL GaussPointsInit
!     thread = 1
! !$    thread = omp_get_thread_num()+1
!     p => IntegStuff(thread)
      p => IntegStuff
      p % n = 1
      p % u(1) = 0
      p % v(1) = 0
      p % w(1) = 0
      p % s(1) = 1
!------------------------------------------------------------------------------
   END FUNCTION GaussPoints0D
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 1D line element
!------------------------------------------------------------------------------
   FUNCTION GaussPoints1D( n ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: n   !< number of points in the requested rule
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
!     INTEGER :: thread, omp_get_thread_num

      IF ( .NOT. GInit ) CALL GaussPointsInit
!     thread = 1
! !$    thread = omp_get_thread_num()+1
!      p => IntegStuff(thread)
      p => IntegStuff
      IF ( n < 1 .OR. n > MAXN ) THEN
        p % n = 0
        WRITE( Message, * ) 'Invalid number of points: ',n
        CALL Error( 'GaussPoints1D', Message )
        RETURN
      END IF

      p % n = n
      p % u(1:n) = Points(1:n,n)
      p % v(1:n) = 0.0d0
      p % w(1:n) = 0.0d0
      p % s(1:n) = Weights(1:n,n)
!------------------------------------------------------------------------------
   END FUNCTION GaussPoints1D
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
   FUNCTION GaussPointsPTriangle(n) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: i,n
      TYPE(GaussIntegrationPoints_t), POINTER :: p
      REAL (KIND=dp) :: uq, vq, sq
!     INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
      IF ( .NOT. GInit ) CALL GaussPointsInit
!      thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      ! Construct Gauss points for p (barycentric) triangle from
      ! Gauss points for quadrilateral
      p = GaussPointsQuad( n )

      ! For each point apply mapping from quad to triangle and
      ! multiply weight by detJ of mapping
!DIR$ IVDEP
      DO i=1,p % n
         uq = p % u(i)
         vq = p % v(i)
         sq = p % s(i)
         p % u(i) = 1d0/2*(uq-uq*vq)
         p % v(i) = SQRT(3d0)/2*(1d0+vq)
         p % s(i) = -SQRT(3d0)/4*(-1+vq)*sq
      END DO

      p % w(1:n) = 0.0d0
!------------------------------------------------------------------------------
    END FUNCTION GaussPointsPTriangle
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 2D triangle element. Here the
!>    reference element over which the integration is done can also be the
!>    equilateral triangle used in the description of p-elements. In that case,
!>    this routine may return a more economical set of integration points.
!------------------------------------------------------------------------------
   FUNCTION GaussPointsTriangle( n, PReferenceElement ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: n    !< number of points in the requested rule
      LOGICAL, OPTIONAL ::  PReferenceElement !< used for switching the reference element
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      INTEGER :: i
      LOGICAL :: ConvertToPTriangle
      REAL (KIND=dp) :: uq, vq, sq
!     INTEGER :: thread, omp_get_thread_num
!-------------------------------------------------------------------------------
      ConvertToPtriangle = .FALSE.
      IF ( PRESENT(PReferenceElement) ) THEN
         ConvertToPTriangle =  PReferenceElement
      END IF

      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      SELECT CASE (n)
      CASE (1)
         p % u(1:n) = UTriangle1P
         p % v(1:n) = VTriangle1P
         p % s(1:n) = STriangle1P / 2.0D0
         p % n = 1
      CASE (3)
         p % u(1:n) = UTriangle3P
         p % v(1:n) = VTriangle3P
         p % s(1:n) = STriangle3P / 2.0D0
         p % n = 3
      CASE (4)
         p % u(1:n) = UTriangle4P
         p % v(1:n) = VTriangle4P
         p % s(1:n) = STriangle4P / 2.0D0
         p % n = 4
      CASE (6)
         p % u(1:n) = UTriangle6P
         p % v(1:n) = VTriangle6P
         p % s(1:n) = STriangle6P / 2.0D0
         p % n = 6
      CASE (7)
         p % u(1:n) = UTriangle7P
         p % v(1:n) = VTriangle7P
         p % s(1:n) = STriangle7P / 2.0D0
         p % n = 7
      CASE (11)
         p % u(1:n) = UTriangle11P
         p % v(1:n) = VTriangle11P
         p % s(1:n) = STriangle11P
         p % n = 11
      CASE (12)
         p % u(1:n) = UTriangle12P
         p % v(1:n) = VTriangle12P
         p % s(1:n) = STriangle12P
         p % n = 12
      CASE (17)
         p % u(1:n) = UTriangle17P
         p % v(1:n) = VTriangle17P
         p % s(1:n) = STriangle17P
         p % n = 17
      CASE (20)
         p % u(1:n) = UTriangle20P
         p % v(1:n) = VTriangle20P
         p % s(1:n) = STriangle20P
         p % n = 20
      CASE DEFAULT
!        CALL Error( 'GaussPointsTriangle', 'Invalid number of points requested' )
!        p % n = 0

        p = GaussPointsQuad( n )
!DIR$ IVDEP
         DO i=1,p % n
            p % v(i) = (p % v(i) + 1) / 2
            p % u(i) = (p % u(i) + 1) / 2 * (1 - p % v(i))
            p % s(i) = p % s(i) * ( 1 - p % v(i) )
         END DO
         p % s(1:p % n) = 0.5d0 * p % s(1:p % n) / SUM( p % s(1:p % n) )
      END SELECT

      IF (ConvertToPTriangle) THEN
         !-------------------------------------------------------------------
         ! Apply an additional transformation if the actual reference element
         ! is the equilateral triangle used in the description of p-elements.
	 ! We map the original integration points into their counterparts on the
	 ! p-reference element and scale the weights by the determinant of the
	 ! deformation gradient associated with the change of reference element.
         !-------------------------------------------------------------------
!DIR$ IVDEP
        DO i=1,P % n
            uq = P % u(i)
            vq = P % v(i)
            sq = P % s(i)
            P % u(i) = -1.0d0 + 2.0d0*uq + vq
            P % v(i) = SQRT(3.0d0)*vq
            P % s(i) = SQRT(3.0d0)*2.0d0*sq
         END DO
      END IF

      p % w(1:n) = 0.0d0
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsTriangle
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!> Return Gaussian integration points for 2D quad element
!------------------------------------------------------------------------------
   FUNCTION GaussPointsQuad( np, PMethod) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: np     !< number of points in the requested rule
      LOGICAL, OPTIONAL ::  PMethod
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      LOGICAL :: Economic
      INTEGER i,j,n,t
!      INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
      Economic = .FALSE.
      IF (PRESENT(PMethod)) Economic = PMethod

      IF ( .NOT. GInit ) CALL GaussPointsInit
!      thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      IF (Economic .AND. (np > 4) .AND. (np <= 60)) THEN
        !PRINT *, 'SELECTING A SPECIAL QUADRATURE FOR p-ELEMENTS'
        SELECT CASE(np)
        CASE (8)
          p % u(1:np) = UPQuad8(1:np)
          p % v(1:np) = VPQuad8(1:np)
          p % s(1:np) = SPQuad8(1:np)
          p % n = 8
        CASE(12)
          p % u(1:np) = UPQuad12(1:np)
          p % v(1:np) = VPQuad12(1:np)
          p % s(1:np) = SPQuad12(1:np)
          p % n = 12
        CASE(20)
          p % u(1:np) = UPQuad20(1:np)
          p % v(1:np) = VPQuad20(1:np)
          p % s(1:np) = SPQuad20(1:np)
          p % n = 20
        CASE(25)
          p % u(1:np) = UPQuad25(1:np)
          p % v(1:np) = VPQuad25(1:np)
          p % s(1:np) = SPQuad25(1:np)
          p % n = 25
        CASE(36)
          p % u(1:np) = UPQuad36(1:np)
          p % v(1:np) = VPQuad36(1:np)
          p % s(1:np) = SPQuad36(1:np)
          p % n = 36
        CASE(45)
          p % u(1:np) = UPQuad45(1:np)
          p % v(1:np) = VPQuad45(1:np)
          p % s(1:np) = SPQuad45(1:np)
          p % n = 45
        CASE(60)
          p % u(1:np) = UPQuad60(1:np)
          p % v(1:np) = VPQuad60(1:np)
          p % s(1:np) = SPQuad60(1:np)
          p % n = 60
        CASE DEFAULT
          WRITE( Message, * ) 'Invalid number of points for p-quadrature: ', np
          CALL Fatal('GaussPointsQuad', Message )
        END SELECT

        p % w(1:np) = 0.0_dp
        !PRINT *, 'NUMBER OF G-POINTS=',p % n
        RETURN
      END IF

      n = SQRT( REAL(np) ) + 0.5

      ! WRITE (*,*) 'Integration:', n, np

      IF ( n < 1 .OR. n > MAXN ) THEN
        p % n = 0
        WRITE( Message, * ) 'Invalid number of points: ', n
        CALL Error( 'GaussPointsQuad', Message )
        RETURN
      END IF

      t = 0
      DO i=1,n
!DIR$ IVDEP
        DO j=1,n
          t = t + 1
          p % u(t) = Points(j,n)
          p % v(t) = Points(i,n)
          p % s(t) = Weights(i,n)*Weights(j,n)
        END DO
      END DO
      p % n = t
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsQuad
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
   FUNCTION GaussPointsPTetra(np) RESULT(p)
!------------------------------------------------------------------------------
   INTEGER :: i,np,n
   TYPE(GaussIntegrationPoints_t), POINTER :: p
   REAL(KIND=dp) :: uh, vh, wh, sh
!  INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
   IF ( .NOT. GInit ) CALL GaussPointsInit
!    thread = 1
! !$ thread = omp_get_thread_num()+1
!    p => IntegStuff(thread)
   p => IntegStuff
   n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0

   ! Get Gauss points of p brick
   ! (take into account term z^2) from jacobian determinant
   p = GaussPointsPBrick(n,n,n+1)
   ! p = GaussPointsBrick( np )
   ! WRITE (*,*) 'Getting Gauss points for: ', n, p % n

   ! For each point apply mapping from brick to
   ! tetrahedron and multiply each weight by detJ
   ! of mapping
!DIR$ IVDEP
   DO i=1,p % n
      uh = p % u(i)
      vh = p % v(i)
      wh = p % w(i)
      sh = p % s(i)

      p % u(i)= 1d0/4*(uh - uh*vh - uh*wh + uh*vh*wh)
      p % v(i)= SQRT(3d0)/4*(5d0/3 + vh - wh/3 - vh*wh)
      p % w(i)= SQRT(6d0)/3*(1d0 + wh)
      p % s(i)= -sh * SQRT(2d0)/16 * (1d0 - vh - wh + vh*wh) * (-1d0 + wh)
   END DO
!------------------------------------------------------------------------------
 END FUNCTION GaussPointsPTetra
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D tetrahedral element. Here the
!>    reference element over which the integration is done can also be the
!>    regular tetrahedron used in the description of p-elements. In that case,
!>    this routine may return a more economical set of integration points.
!------------------------------------------------------------------------------
   FUNCTION GaussPointsTetra( n, PReferenceElement ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: n      !< number of points in the requested rule
      LOGICAL, OPTIONAL ::  PReferenceElement !< used for switching the reference element
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      REAL( KIND=dp ) :: ScaleFactor
      INTEGER :: i
      LOGICAL :: ConvertToPTetrahedron
      REAL (KIND=dp) :: uq, vq, wq, sq
!     INTEGER :: thread, omp_get_thread_num
!----------------------------------------------------------------------------------
      ConvertToPTetrahedron = .FALSE.
      IF ( PRESENT(PReferenceElement) ) THEN
         ConvertToPTetrahedron =  PReferenceElement
      END IF

      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      SELECT CASE (n)
      CASE (1)
         p % u(1:n) = UTetra1P
         p % v(1:n) = VTetra1P
         p % w(1:n) = WTetra1P
         p % s(1:n) = STetra1P / 6.0D0
         p % n = 1
      CASE (4)
         p % u(1:n) = UTetra4P
         p % v(1:n) = VTetra4P
         p % w(1:n) = WTetra4P
         p % s(1:n) = STetra4P / 6.0D0
         p % n = 4
      CASE (5)
         p % u(1:n) = UTetra5P
         p % v(1:n) = VTetra5P
         p % w(1:n) = WTetra5P
         p % s(1:n) = STetra5P / 6.0D0
         p % n = 5
      CASE (11)
         p % u(1:n) = UTetra11P
         p % v(1:n) = VTetra11P
         p % w(1:n) = WTetra11P
         p % s(1:n) = STetra11P / 6.0D0
         p % n = 11
      CASE DEFAULT
!        CALL Error( 'GaussPointsTetra', 'Invalid number of points requested.' )
!        p % n = 0
         p = GaussPointsBrick( n )

!DIR$ IVDEP
         DO i=1,p % n
            ScaleFactor = 0.5d0
            p % u(i) = ( p % u(i) + 1 ) * Scalefactor
            p % v(i) = ( p % v(i) + 1 ) * ScaleFactor
            p % w(i) = ( p % w(i) + 1 ) * ScaleFactor
            p % s(i) = p % s(i) * ScaleFactor**3

            ScaleFactor = 1.0d0 - p % w(i)
            p % u(i) = p % u(i) * ScaleFactor
            p % v(i) = p % v(i) * ScaleFactor
            p % s(i) = p % s(i) * ScaleFactor**2

            ScaleFactor = 1.0d0 - p % v(i) / ScaleFactor
            p % u(i) = p % u(i) * ScaleFactor
            p % s(i) = p % s(i) * ScaleFactor
         END DO
!         p % s(1:p % n) = p % s(1:p % n) / SUM( p % s(1:p % n) ) / 6.0d0
      END SELECT

      IF (ConvertToPTetrahedron) THEN
         !-------------------------------------------------------------------
         ! Apply an additional transformation if the actual reference element
         ! is the regular tetrahedron used in the description of p-elements
	 ! We map the original integration points into their counterparts on the
	 ! p-reference element and scale the weights by the determinant of the
	 ! deformation gradient associated with the change of reference element.
         !-------------------------------------------------------------------
!DIR$ IVDEP
        DO i=1,P % n
            uq = P % u(i)
            vq = P % v(i)
            wq = P % w(i)
            sq = P % s(i)
            P % u(i) = -1.0d0 + 2.0d0*uq + vq + wq
            P % v(i) = SQRT(3.0d0)*vq + 1.0d0/SQRT(3.0d0)*wq
            P % w(i) = SQRT(8.0d0)/SQRT(3.0d0)*wq
            P % s(i) = SQRT(8.0d0)*2.0d0*sq
         END DO
      END IF
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsTetra
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
   FUNCTION GaussPointsPPyramid( np ) RESULT(p)
!------------------------------------------------------------------------------
   INTEGER :: np,n,i
   REAL(KIND=dp) :: uh,vh,wh,sh
   TYPE(GaussIntegrationPoints_t), POINTER :: p
!  INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
   IF ( .NOT. GInit ) CALL GaussPointsInit
!    thread = 1
! !$ thread = omp_get_thread_num()+1
!    p => IntegStuff(thread)
   p => IntegStuff

   n = DBLE(np)**(1.0D0/3.0D0) + 0.5D0

   ! Get Gauss points of p brick
   ! (take into account term (-1+z)^2) from jacobian determinant
   p = GaussPointsPBrick(n,n,n+1)

   ! For each point apply mapping from brick to
   ! pyramid and multiply each weight by detJ
   ! of mapping
!DIR$ IVDEP
   DO i=1,p % n
      uh = p % u(i)
      vh = p % v(i)
      wh = p % w(i)
      sh = p % s(i)

      p % u(i)= 1d0/2*uh*(1d0-wh)
      p % v(i)= 1d0/2*vh*(1d0-wh)
      p % w(i)= SQRT(2d0)/2*(1d0+wh)
      p % s(i)= sh * SQRT(2d0)/8 * (-1d0+wh)**2
   END DO
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsPPyramid
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D pyramid element.
!------------------------------------------------------------------------------
   FUNCTION GaussPointsPyramid( np ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: np     !< number of points in the requested rule
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      INTEGER :: i,j,k,n,t
!       INTEGER :: thread, omp_get_thread_num

      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      n = REAL(np)**(1.0D0/3.0D0) + 0.5D0

      IF ( n < 1 .OR. n > MAXN ) THEN
         p % n = 0
         WRITE( Message, * ) 'Invalid number of points: ', n
         CALL Error( 'GaussPointsPyramid', Message )
         RETURN
      END IF

      t = 0
      DO i=1,n
        DO j=1,n
!DIR$ IVDEP
          DO k=1,n
             t = t + 1
             p % u(t) = Points(k,n)
             p % v(t) = Points(j,n)
             p % w(t) = Points(i,n)
             p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
          END DO
        END DO
      END DO
      p % n = t

!DIR$ IVDEP
      DO t=1,p % n
        p % w(t) = (p % w(t) + 1.0d0) / 2.0d0
        p % u(t) = p % u(t) * (1.0d0-p % w(t))
        p % v(t) = p % v(t) * (1.0d0-p % w(t))
        p % s(t) = p % s(t) * (1.0d0-p % w(t))**2/2
      END DO
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsPyramid
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
   FUNCTION GaussPointsPWedge(n) RESULT(p)
!------------------------------------------------------------------------------
   INTEGER :: n, i
   REAL(KIND=dp) :: uh,vh,wh,sh
   TYPE(GaussIntegrationPoints_t), POINTER :: p
!   INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
   IF ( .NOT. GInit ) CALL GaussPointsInit
!    thread = 1
! !$ thread = omp_get_thread_num()+1
!    p => IntegStuff(thread)
   p => IntegStuff

   ! Get Gauss points of brick
   p = GaussPointsBrick(n)

   ! For each point apply mapping from brick to
   ! wedge and multiply each weight by detJ
   ! of mapping
!DIR$ IVDEP
   DO i=1,p % n
      uh = p % u(i)
      vh = p % v(i)
      wh = p % w(i)
      sh = p % s(i)

      p % u(i)= 1d0/2*(uh-uh*vh)
      p % v(i)= SQRT(3d0)/2*(1d0+vh)
      p % w(i)= wh
      p % s(i)= sh * SQRT(3d0)*(1-vh)/4
   END DO
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsPWedge
!------------------------------------------------------------------------------


   
!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D wedge element
!------------------------------------------------------------------------------
   FUNCTION GaussPointsWedge( np ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: np     !< number of points in the requested rule
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      INTEGER :: i,j,k,n,t
!       INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      n = REAL(np)**(1.0d0/3.0d0) + 0.5d0

      IF ( n < 1 .OR. n > MAXN ) THEN
         p % n = 0
         WRITE( Message, * ) 'Invalid number of points: ', n
         CALL Error( 'GaussPointsWedge', Message )
         RETURN
      END IF

      t = 0
      DO i=1,n
        DO j=1,n
!DIR$ IVDEP
          DO k=1,n
             t = t + 1
             p % u(t) = Points(k,n)
             p % v(t) = Points(j,n)
             p % w(t) = Points(i,n)
             p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
          END DO
        END DO
      END DO
      p % n = t

!DIR$ IVDEP
      DO i=1,p % n
        p % v(i) = (p % v(i) + 1)/2
        p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i))
        p % s(i) = p % s(i) * (1-p % v(i))/4
      END DO
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsWedge
!------------------------------------------------------------------------------

!------------------------------------------------------------------------------
!>  Return an optimized number of Gaussian points for integrating over prisms.
!>  Here the reference element can also be that of the p-approximation.
!>  A rule with m x n points is returned.
!------------------------------------------------------------------------------
   FUNCTION GaussPointsWedge2(m,n,PReferenceElement) RESULT(p)
!------------------------------------------------------------------------------
      TYPE(GaussIntegrationPoints_t), POINTER :: p
      INTEGER :: m     !< The number of points over a triangular face
      INTEGER :: n     !< The number of points in the orthogonal direction to the triangular faces
      LOGICAL, OPTIONAL ::  PReferenceElement !< Used for switching the reference element
!------------------------------------------------------------------------------
      INTEGER :: i,j,k,t
      LOGICAL :: ConvertToPPrism
      REAL(KIND=dp) :: uq(20), vq(20), sq(20)
      REAL(KIND=dp) :: uh,vh,wh,sh
!-----------------------------------------------------------------------------
      ConvertToPPrism = .FALSE.
      IF ( PRESENT(PReferenceElement) ) THEN
         ConvertToPPrism =  PReferenceElement
      END IF
      IF ( .NOT. GInit ) CALL GaussPointsInit
      p => IntegStuff

      SELECT CASE (m)
      CASE (1)
         IF (ConvertToPPrism) THEN
            uq(1) = -1.0d0 + 2.0d0*UTriangle1P(1) + VTriangle1P(1)
            vq(1) = SQRT(3.0d0) * VTriangle1P(1)
            sq(1) = SQRT(3.0d0) * STriangle1P(1)
         ELSE
            uq(1:m) = UTriangle1P
            vq(1:m) = VTriangle1P
            sq(1:m) = STriangle1P / 2.0D0
         END IF
      CASE (3)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle3P(i) + VTriangle3P(i)
               vq(i) = SQRT(3.0d0) * VTriangle3P(i)
               sq(i) = SQRT(3.0d0) * STriangle3P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle3P
            vq(1:m) = VTriangle3P
            sq(1:m) = STriangle3P / 2.0D0
         END IF
      CASE (4)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle4P(i) + VTriangle4P(i)
               vq(i) = SQRT(3.0d0) * VTriangle4P(i)
               sq(i) = SQRT(3.0d0) * STriangle4P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle4P
            vq(1:m) = VTriangle4P
            sq(1:m) = STriangle4P / 2.0D0
         END IF
      CASE (6)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle6P(i) + VTriangle6P(i)
               vq(i) = SQRT(3.0d0) * VTriangle6P(i)
               sq(i) = SQRT(3.0d0) * STriangle6P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle6P
            vq(1:m) = VTriangle6P
            sq(1:m) = STriangle6P / 2.0D0
         END IF
      CASE (7)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle7P(i) + VTriangle7P(i)
               vq(i) = SQRT(3.0d0) * VTriangle7P(i)
               sq(i) = SQRT(3.0d0) * STriangle7P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle7P
            vq(1:m) = VTriangle7P
            sq(1:m) = STriangle7P / 2.0D0
         END IF
      CASE (11)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle11P(i) + VTriangle11P(i)
               vq(i) = SQRT(3.0d0) * VTriangle11P(i)
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle11P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle11P
            vq(1:m) = VTriangle11P
            sq(1:m) = STriangle11P
         END IF
      CASE (12)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle12P(i) + VTriangle12P(i)
               vq(i) = SQRT(3.0d0) * VTriangle12P(i)
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle12P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle12P
            vq(1:m) = VTriangle12P
            sq(1:m) = STriangle12P
         END IF
      CASE (17)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle17P(i) + VTriangle17P(i)
               vq(i) = SQRT(3.0d0) * VTriangle17P(i)
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle17P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle17P
            vq(1:m) = VTriangle17P
            sq(1:m) = STriangle17P
         END IF
      CASE (20)
        IF (ConvertToPPrism) THEN
!DIR$ IVDEP
            DO i=1,m
               uq(i) = -1.0d0 + 2.0d0*UTriangle20P(i) + VTriangle20P(i)
               vq(i) = SQRT(3.0d0) * VTriangle20P(i)
               sq(i) = SQRT(3.0d0) * 2.0d0 * STriangle20P(i)
            END DO
         ELSE
            uq(1:m) = UTriangle20P
            vq(1:m) = VTriangle20P
            sq(1:m) = STriangle20P
         END IF
      CASE DEFAULT
         !-------------------------------------------------------------------------------------
         ! The requested number of integration points associated with the triangular face
         ! cannot be created. Therefore use the given number of points in the third direction
         ! and generate silently n x n x n rule.
         !-----------------------------------------------------------------------------------
         IF (ConvertToPPrism) THEN
           p = GaussPointsBrick(n*n*n)
!DIR$ IVDEP
            DO i=1,p % n
               uh = p % u(i)
               vh = p % v(i)
               wh = p % w(i)
               sh = p % s(i)

               p % u(i)= 1d0/2*(uh-uh*vh)
               p % v(i)= SQRT(3d0)/2*(1d0+vh)
               p % w(i)= wh
               p % s(i)= sh * SQRT(3d0)*(1-vh)/4
            END DO
         ELSE
            t = 0
            DO i=1,n
              DO j=1,n
!DIR$ IVDEP
                  DO k=1,n
                     t = t + 1
                     p % u(t) = Points(k,n)
                     p % v(t) = Points(j,n)
                     p % w(t) = Points(i,n)
                     p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
                  END DO
               END DO
            END DO
            p % n = t

!DIR$ IVDEP
            DO i=1,p % n
               p % v(i) = (p % v(i) + 1)/2
               p % u(i) = (p % u(i) + 1)/2 * (1 - p % v(i))
               p % s(i) = p % s(i) * (1-p % v(i))/4
            END DO
         END IF
         RETURN
      END SELECT

      t = 0
      DO i=1,m
!DIR$ IVDEP
         DO j=1,n
            t = t + 1
            p % u(t) = uq(i)
            p % v(t) = vq(i)
            p % w(t) = Points(j,n)
            p % s(t) = sq(i)*Weights(j,n)
         END DO
      END DO
      p % n = t
!------------------------------------------------------------------------------
    END FUNCTION GaussPointsWedge2
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D wedge elements using
!> economical quadratures that are not product of segment and triangle rules.
!------------------------------------------------------------------------------
   FUNCTION GaussPointsWedgeEconomic( n, PReferenceElement ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: n      !< number of points in the requested rule
      LOGICAL, OPTIONAL ::  PReferenceElement !< used for switching the reference element 
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      REAL( KIND=dp ) :: ScaleFactor
      INTEGER :: i
      LOGICAL :: ConvertToPWedge
      REAL (KIND=dp) :: uq, vq, wq, sq
!     INTEGER :: thread, omp_get_thread_num
!----------------------------------------------------------------------------------
      ConvertToPWedge = .FALSE.
      IF ( PRESENT(PReferenceElement) ) THEN
        ConvertToPWedge = PReferenceElement
      END IF

      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      SELECT CASE (n)
      CASE (4)
         p % u(1:n) = UWedge4P
         p % v(1:n) = VWedge4P
         p % w(1:n) = WWedge4P
         p % s(1:n) = SWedge4P
      CASE (5)
         p % u(1:n) = UWedge5P
         p % v(1:n) = VWedge5P
         p % w(1:n) = WWedge5P
         p % s(1:n) = SWedge5P 
      CASE (6)
         p % u(1:n) = UWedge6P
         p % v(1:n) = VWedge6P
         p % w(1:n) = WWedge6P
         p % s(1:n) = SWedge6P 
      CASE (7)
         p % u(1:n) = UWedge7P
         p % v(1:n) = VWedge7P
         p % w(1:n) = WWedge7P
         p % s(1:n) = SWedge7P
      CASE (10)
         p % u(1:n) = UWedge10P
         p % v(1:n) = VWedge10P
         p % w(1:n) = WWedge10P
         p % s(1:n) = SWedge10P 
      CASE (11)
         p % u(1:n) = UWedge11P
         p % v(1:n) = VWedge11P
         p % w(1:n) = WWedge11P
         p % s(1:n) = SWedge11P 
      CASE (14)
         p % u(1:n) = UWedge14P
         p % v(1:n) = VWedge14P
         p % w(1:n) = WWedge14P
         p % s(1:n) = SWedge14P 
      CASE (15)
         p % u(1:n) = UWedge15P
         p % v(1:n) = VWedge15P
         p % w(1:n) = WWedge15P
         p % s(1:n) = SWedge15P 
      CASE (16)
         p % u(1:n) = UWedge16P
         p % v(1:n) = VWedge16P
         p % w(1:n) = WWedge16P
         p % s(1:n) = SWedge16P
      CASE (24)
         p % u(1:n) = UWedge24P
         p % v(1:n) = VWedge24P
         p % w(1:n) = WWedge24P
         p % s(1:n) = SWedge24P
         
      CASE DEFAULT
        CALL Fatal( 'GaussPointsWedgeEconomic',&
            'Invalid number of points requested')
      END SELECT

      p % n = n

      IF (ConvertToPWedge ) THEN
        DO i=1,P % n  
          uq = P % u(i) 
          vq = P % v(i)
          sq = P % s(i)

          ! first to classical
          uq = (uq+1.d0)/2.0d0
          vq = (vq+1.d0)/2.0d0

          ! then to p-convention
          P % u(i) = -1.0d0 + 2.0d0*uq + vq
          P % v(i) = SQRT(3.0d0) * vq            
          P % s(i) = SQRT(3.0d0) * sq / 2.0d0
        END DO
      ELSE
        ! Map to classical Elmer local coordinates in [0,1]
        p % u(1:n) = ( p % u(1:n)+1.0d0 ) / 2.0d0
        p % v(1:n) = ( p % v(1:n)+1.0d0 ) / 2.0d0 
        p % s(1:n) = p % s(1:n) / 4.0d0
      END IF
!------------------------------------------------------------------------------
    END FUNCTION GaussPointsWedgeEconomic
!------------------------------------------------------------------------------

    

!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D brick element for
!>    composite rule
!>    sum_i=1^nx(sum_j=1^ny(sum_k=1^nz w_ijk f(x_{i,j,k},y_{i,j,k},z_{i,j,k}))).
!------------------------------------------------------------------------------
   FUNCTION GaussPointsPBrick( nx, ny, nz ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: nx    !< number of points in the requested rule in x direction
      INTEGER :: ny    !< number of points in the requested rule in y direction
      INTEGER :: nz    !< number of points in the requested rule in z direction
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      INTEGER i,j,k,t
!       INTEGER :: thread, omp_get_thread_num
!------------------------------------------------------------------------------
      IF ( .NOT. GInit ) CALL GaussPointsInit
!       thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      ! Check validity of number of integration points
      IF ( nx < 1 .OR. nx > MAXN .OR. &
           ny < 1 .OR. ny > MAXN .OR. &
           nz < 1 .OR. nz > MAXN) THEN
        p % n = 0
        WRITE( Message, * ) 'Invalid number of points: ', nx, ny, nz
        CALL Error( 'GaussPointsBrick', Message )
        RETURN
      END IF

      t = 0
      DO i=1,nx
        DO j=1,ny
!DIR$ IVDEP
          DO k=1,nz
            t = t + 1
            p % u(t) = Points(i,nx)
            p % v(t) = Points(j,ny)
            p % w(t) = Points(k,nz)
            p % s(t) = Weights(i,nx)*Weights(j,ny)*Weights(k,nz)
          END DO
        END DO
      END DO
      p % n = t
!------------------------------------------------------------------------------
    END FUNCTION GaussPointsPBrick
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>    Return Gaussian integration points for 3D brick element
!------------------------------------------------------------------------------
   FUNCTION GaussPointsBrick( np ) RESULT(p)
!------------------------------------------------------------------------------
      INTEGER :: np     !< number of points in the requested rule
      TYPE(GaussIntegrationPoints_t), POINTER :: p
!------------------------------------------------------------------------------
      INTEGER i,j,k,n,t
!      INTEGER :: thread, omp_get_thread_num

      IF ( .NOT. GInit ) CALL GaussPointsInit
!      thread = 1
! !$    thread = omp_get_thread_num()+1
!       p => IntegStuff(thread)
      p => IntegStuff

      n = REAL(np)**(1.0D0/3.0D0) + 0.5D0

      IF ( n < 1 .OR. n > MAXN ) THEN
        p % n = 0
        WRITE( Message, * ) 'Invalid number of points: ', n
        CALL Error( 'GaussPointsBrick', Message )
        RETURN
      END IF

      t = 0
      DO i=1,n
        DO j=1,n
!DIR$ IVDEP
          DO k=1,n
            t = t + 1
            p % u(t) = Points(k,n)
            p % v(t) = Points(j,n)
            p % w(t) = Points(i,n)
            p % s(t) = Weights(i,n)*Weights(j,n)*Weights(k,n)
          END DO
        END DO
      END DO
      p % n = t
!------------------------------------------------------------------------------
   END FUNCTION GaussPointsBrick
!------------------------------------------------------------------------------


!------------------------------------------------------------------------------
!>    Given element structure return Gauss integration points for the element.
!----------------------------------------------------------------------------------------------
   FUNCTION GaussPoints( elm, np, RelOrder, EdgeBasis, PReferenceElement, &
        EdgeBasisDegree) RESULT(IntegStuff)
!---------------------------------------------------------------------------------------------
     USE PElementMaps, ONLY : isActivePElement
     TYPE( Element_t ) :: elm   !< Structure holding element description
     INTEGER, OPTIONAL :: np    !< Number of requested points
     INTEGER, OPTIONAL :: RelOrder  !< Relative order of Gauss points (-1,0,+1)
     LOGICAL, OPTIONAL :: EdgeBasis !< For choosing quadrature suitable for edge elements
     LOGICAL, OPTIONAL :: PReferenceElement !< For switching to the p-version reference element
     INTEGER, OPTIONAL :: EdgeBasisDegree ! The degree of edge elements
     TYPE( GaussIntegrationPoints_t ) :: IntegStuff   !< Structure holding the integration points
!------------------------------------------------------------------------------
     LOGICAL :: pElement, UsePRefElement, Economic
     INTEGER :: n, eldim, p1d, ntri, nseg, necon
     TYPE(ElementType_t), POINTER :: elmt
!------------------------------------------------------------------------------
     elmt => elm % TYPE

     IF (PRESENT(EdgeBasis)) THEN
       IF (EdgeBasis) THEN
         UsePRefElement = .TRUE.
         IF (PRESENT(PReferenceElement)) UsePRefElement = PReferenceElement
         IF (PRESENT(EdgeBasisDegree)) THEN
           IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement, EdgeBasisDegree)
         ELSE
           IntegStuff = EdgeElementGaussPoints(elmt % ElementCode/100, UsePRefElement)
         END IF
         RETURN
       END IF
     END IF

     IF( PRESENT(PReferenceElement)) THEN
       pElement = PReferenceElement
     ELSE
       pElement = isActivePElement(elm)
     END IF
     
     IF ( PRESENT(np) ) THEN
       n = np
     ELSE IF( PRESENT( RelOrder ) ) THEN
       IF (pElement) THEN
         n = elm % PDefs % GaussPoints
         IF( RelOrder == 0 ) THEN
           CONTINUE
         ELSE
           eldim = elm % type % dimension
           p1d = NINT( n**(1.0_dp/eldim) ) + RelOrder
           IF( p1d < 1 ) THEN
             CALL Fatal('GaussPoints','Number of integration points must remain positive!')
           END IF
           n = p1d**eldim

           ! Take into account the case of economic integration. Map a tensor
           ! product rule to a corresponding economic rule.
           IF (elmt % ElementCode / 100 == 4) THEN
             SELECT CASE(n)
             CASE(9)
               n = 8
             CASE(16)
               n = 12
             CASE(25)
               n = 20
             CASE(36)
               n = 25
             CASE(49)
               n = 36
             CASE(64)
               n = 45
             CASE(81)
               n = 60
             CASE DEFAULT
               CONTINUE
             END SELECT
           END IF
         END IF
       ELSE
         IF( RelOrder == 0 ) THEN
           n = elmt % GaussPoints
         ELSE IF( RelOrder == 1 ) THEN
           n = elmt % GaussPoints2
         ELSE IF( RelOrder == -1 ) THEN
           n = elmt % GaussPoints0
         ELSE
           PRINT *,'RelOrder can only be {-1, 0, 1} !'
         END IF
       END IF
     ELSE
       IF (pElement) THEN
         IF( ASSOCIATED( elm % PDefs ) ) THEN
           n = elm % PDefs % GaussPoints
         ELSE
           n = elmt % GaussPoints
         END IF
         IF( n == 0 ) n = elmt % GaussPoints
       ELSE
         n = elmt % GaussPoints
       END IF
     END IF

     SELECT CASE( elmt % ElementCode / 100 )
     CASE (1)
        IntegStuff = GaussPoints0D(n)

     CASE (2)
        IntegStuff = GaussPoints1D(n)

     CASE (3)
        IF (pElement) THEN
          IntegStuff = GaussPointsPTriangle(n)
        ELSE
          IntegStuff = GaussPointsTriangle(n)
        END IF

     CASE (4)
       IF (pElement .AND. ASSOCIATED( elm % pdefs ) ) THEN
         Economic = .FALSE.
         ! For certain polynomial orders, economic quadratures may be used:
         IF (elm % PDefs % P > 1 .AND.  elm % PDefs % P <= 8) Economic = .TRUE.
         ! An explicit bubble augmentation with lower-order methods switches to
         ! the standard rule:
         IF (elm % BDOFs > 0 .AND. elm % PDefs % P < 4) Economic = .FALSE.
         IntegStuff = GaussPointsQuad(n, Economic)
       ELSE
         IntegStuff = GaussPointsQuad(n)
       END IF

     CASE (5)
        IF (pElement) THEN
           IntegStuff = GaussPointsPTetra(n)
        ELSE
           IntegStuff = GaussPointsTetra(n)
        END IF

     CASE (6)
        IF (pElement) THEN
           IntegStuff = GaussPointsPPyramid(n)
        ELSE
           IntegStuff = GaussPointsPyramid(n)
        END IF

      CASE (7)
        IF( PRESENT( np ) ) THEN
          ! possible values:
          ! triangle = 1, 3, 4, 6, 7, 11, 12, 17, 20
          ! segment  = 1, 2, 3, 4, 5, 6,  7,  8,  9,

          ntri = 0; nseg = 0; necon = 0

          SELECT CASE( n )

            ! Cases ordered so that we have the segment x triangle rules first
          CASE( 1, 2, 3)
            nseg = n
          CASE( 6, 8 )
            nseg = 2
          CASE( 12, 18, 21 )
            nseg = 3
          CASE( 28, 44, 48 )
            nseg = 4
          CASE( 85, 100 )
            nseg = 5

            ! The economical rules
          CASE( 4, 5, 7, 10, 11, 14, 15, 16, 24 )
            ! Note: we would have 6 and 8 point rules from the economic family as well
            necon = n 
          END SELECT

          IF( nseg > 0 ) THEN
            ntri =  n / nseg
            IntegStuff = GaussPointsWedge2(ntri,nseg,PReferenceElement=pElement)
            RETURN
          ELSE IF( necon > 0 ) THEN
            IntegStuff = GaussPointsWedgeEconomic(necon,PReferenceElement=pElement)
            RETURN
          END IF
        END IF

        IF (pElement) THEN
           IntegStuff = GaussPointsPWedge(n)
        ELSE
           IntegStuff = GaussPointsWedge(n)
        END IF

     CASE (8)
        IntegStuff = GaussPointsBrick(n)
     END SELECT

   END FUNCTION GaussPoints


!----------------------------------------------------------------------------------
!>  Return a suitable version of the Gaussian numerical integration method for
!>  H(curl)-conforming finite elements. The default here is that the edge basis
!>  functions are obtained via calling the function EdgeElementInfo, so that
!>  the reference element is chosen to be that used for p-approximation.
!>  By giving the optional argument PiolaVersion = .FALSE. this function returns
!>  an alternate (usually smaller) set of integration points on the classic
!>  reference element. This option may be useful when the edge basis functions are
!>  obtained via calling the alternate subroutine GetEdgeBasis.
!----------------------------------------------------------------------------------------
   FUNCTION EdgeElementGaussPoints(ElementFamily, PiolaVersion, BasisDegree) RESULT(IP)
!---------------------------------------------------------------------------------------
     INTEGER :: ElementFamily
     LOGICAL, OPTIONAL :: PiolaVersion
     INTEGER, OPTIONAL :: BasisDegree
     TYPE(GaussIntegrationPoints_t) :: IP
!------------------------------------------------------------------------------
     LOGICAL :: PRefElement, SecondOrder
!------------------------------------------------------------------------------
     PRefElement = .TRUE.
     SecondOrder = .FALSE.
     IF ( PRESENT(PiolaVersion) ) PRefElement = PiolaVersion
     IF ( PRESENT(BasisDegree) ) SecondOrder = BasisDegree > 1

     SELECT CASE(ElementFamily)
     CASE (1)
        IP = GaussPoints0D(1)
     CASE (2)
        IP = GaussPoints1D(2)
     CASE(3)
        IF (SecondOrder) THEN
           IP = GaussPointsTriangle(6, PReferenceElement=PRefElement)
        ELSE
           IP = GaussPointsTriangle(3, PReferenceElement=PRefElement)
        END IF
     CASE(4)
        IF (PRefElement) THEN
           IP = GaussPointsQuad(9)
        ELSE
           IP = GaussPointsQuad(4)
        END IF
     CASE(5)
        IF (SecondOrder) THEN
           IP = GaussPointsTetra(11, PReferenceElement=PRefElement)
        ELSE
           IP = GaussPointsTetra(4, PReferenceElement=PRefElement)
        END IF
     CASE(6)
        IF (PRefElement) THEN
           IP = GaussPointsPPyramid(27)
        ELSE
           IP = GaussPointsPyramid(27)
        END IF
     CASE(7)
        IF (PRefElement) THEN
           IP = GaussPointsWedge2(6,3,PReferenceElement=PRefElement)
        ELSE
           IP = GaussPointsWedge2(3,2,PReferenceElement=PRefElement)
        END IF
     CASE(8)
        IF (PRefElement) THEN
           IP = GaussPointsPBrick(3,3,3)
        ELSE
           IP = GaussPointsBrick(8)
        END IF
     CASE DEFAULT
        CALL Fatal('Integration::EdgeElementGaussPoints','Unsupported element type')
     END SELECT
!------------------------------------------------------------------------------
   END FUNCTION EdgeElementGaussPoints
!------------------------------------------------------------------------------


!---------------------------------------------------------------------------
END MODULE Integration
!---------------------------------------------------------------------------

!> \} ElmerLib
